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This  report  was  prepared  fay  J.  E.  Edinger  Associates,  Inc.,  (JEEAI) 
for  the  U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES)  under 
Purchase  Order  No.  DACW29-81-M-2788  dated  6  May  1981.  The  study  was 
sponsored  by  the  Office,  Chief  of  Engineers,  under  the  Environmental 
Impact  Research  Program. 

The  study  was  monitored  at  WES  by  Mr.  Ross  W.  Hall,  Jr.,  Ecosystem 
Research  and  Simulation  Division  (ERSD),  Environmental  Laboratory  (EL). 

Mr.  Donald  L.  Robey,  Chief,  ERSD,  and  Dr.  John  Harrison,  Chief,  EL,  pro¬ 
vided  general  supervision. 

Director  of  WES  during  report  preparation  was  Col.  Tilford  C.  Creel, 

CE.  The  Technical  Director  was  Mr.  F.  R.  Brown. 

This  report  should  be  referenced  as  follows: 

Buchak,  E.M.,  and  Edinger,  J.E.  1982.  "User  Guide  for 
CE-QUAL-ELV2 :  A  Longitudinal-Vertical,  Time-Varying 
Estuarine  Water  Quality  Model,"  Instruction  Report 
EL-82-1,  U.  S.  Army  Engineer  Waterways  Experiment  Sta¬ 
tion,  CE,  Vicksburg,  Miss. 
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CONVERSION  FACTORS,  INCH-POUND  TO  METRIC 
(SI)  UNITS  OF  MEASUREMENT 


Inch-pound  units  of  measurement  used  in  this  report  can  be  converted  to 
metric  (SI)  units  as  follows: 


Multiply 

By 

To  Obtain 

inches 

0.0254 

meters 

miles  (international 

nautical) 

1852.0 

meters 

feet 

0.3048 

meters 

Btu  ( thermochemical) / 

( clay  f 2  •  °F) 

0.236436 

watts  per  square  meter  Kelvin 

Fahrenheit  degrees 

5/9 

Celsius  degrees  or  Kelvins* 

*  T„  obtain  Celsius  (C)  temperature  readings  from  Fahrenheit  (F)  readings, 
use  the  following  formula:  C  =  (5/9) (F  -  32).  To  obtain  Kelvin  (K) 
readings,  use:  K  =  (5/9) (F  -  32)  +  273.15  . 


USER  GUIDE  FOR  CE-QUAL-ELV2 :  A  LONGITUDINAL-VERTICAL 


TIME- VARYING  ESTUARINE  WATER  QUALITY  MODEL 
PART  I:  CAPABILITIES  AND  LIMITATIONS 


1.  CE-QUAL-ELV2  was  developed  to  assist  in  the  analysis  of  water  qual¬ 
ity  problems  in  estuaries  where  buoyancy  is  important  and  where  lateral 
homogeneity  can  be  assumed.  The  code  is  recommended  for  those  cases  where 
longitudinal  and  vertical  temperature,  salinity,  or  constituent  gradients 
occur.  CE-QUAL-ELV2  generates  time-varying  velocity,  temperature,  salin¬ 
ity,  and  water  quality  constituent  fields  and  surface  elevations  on  a 
longitudinal  and  vertical  grid.  Both  the  spatial  and  temporal  resolution 
can  be  varied  by  the  user,  within  certain  limits.  Units  of  measure  are 
expressed  according  to  the  International  System  of  Units  (SI) ,  (meter- 
kilogram-second) . 

2.  An  application  requires  a  considerable  effort  in  planning  and  data 
assembly  and  computer  resources.  The  user  must  modify  a  subroutine  that 
supplies  time-varying  boundary  condition  data  to  CE-QUAL-ELV2  for  each  par¬ 
ticular  case.  If  the  user  intends  to  make  use  of  the  constituent  transport 
computation  option,  he  must  be  able  to  quantify  the  reaction  rates  for  each 
constituent  in  terms  of  every  other  constituent  and  code  these  statements. 
The  user  must  also  be  able  to  review  results  for  reasonableness  and  appli¬ 
cability  to  his  problem. 

3.  CE-QUAL-ELV2  in  its  present  form  can  be  applied  to  a  single,  con¬ 
tinuous  reach  of  an  estuary.  Carrying  the  computations  into  branches  and 
tributaries  is  possible  if  the  code  is  employed  as  a  subroutine  for  each 
reach  and  boundary  conditions  between  reaches  are  specified.  This  config¬ 
uration  has  not  been  tested,  however,  and  would  require  significant  addi¬ 
tional  programming. 

4.  The  code  automatically  increases  the  number  of  horizontal  layers 
when  water  surface  elevations  increase  and  then  adds  segments  as  back¬ 
waters  progress  into  the  freshwater  end  of  the  estuary.  Similarly,  the 
grid  contracts  both  vertically  and  longitudinally  when  surface  elevations 
decrease . 


5.  Simulations  of  long  periods  are  economically  feasible.  An  eight- 
month  simulation  at  a  time  step  of  15  minutes  for  temperature  and  salinity 
on  a  grid  that  is  30  segments  by  20  layers  with  375  active  cells  would 
cost  $300  on  a  commerical  batch  processing  computer  system.  Each  water 
quality  constituent  would  add  approximately  $30  to  this  cost.  These  cost 
estimates  are  exclusive  of  any  intermediate  simulations  and  postprocessing 
of  results.  Experience  has  shown  that  ten  intermediate  simulations  may  be 
required  for  each  final  simulation  obtained. 

b.  CE-QUAL-ELV2  has  been  tested  against  analytic  solutions  and  has  been 
applied  to  several  cases.  The  program  from  which  CE-QUAL-ELV2  is  derived, 
the  Laterally  Averaged  Reservoir  Model,  Version  2  (LARM2) ,  has  been  success¬ 
fully  applied  to  approximately  twenty  cases. 
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PART  II:  THEORETICAL  BASIS,  ESTUARINE  BOUNDARY 
CONDITIONS,  AND  VALIDATION  TESTS 


7.  Details  of  the  formulation  of  the  governing  partial  differential 
equations  and  the  subsequent  casting  of  those  equations  into  finite  dif¬ 
ference  form  are  presented  in  Edinger  and  Buchak  (1979  and  1980a)  . 

Seven  equations  (longitudinal  momentum,  vertical  momentum,  continuity, 
heat,  salinity,  and  constituent  balances,  and  state)  are  solved  to  ob¬ 
tain  longitudinal  and  vertical  velocity  components,  surface  elevations, 
temperatures,  salinities,  constituent  concentrations,  and  densities  at  a 
grid  of  points  in  space  and  time.  The  three-dimensional,  time-varying 
partial  differential  equations  are  formally  averaged  over  the  estuary 
width  and  then  cast  into  finite  difference  form  after  vertical  averaging 
over  a  horizontal  layer  thickness  h  with  the  k  ^  layer  boundaries  at 
z  =  k  +  h  and  z  ■  k  -  h  .  The  laterally  averaged  equations  as  verti¬ 
cally  integrated  over  the  layer  thickness  are  presented  below: 


longitudinal  (x-direction)  momentum 


k  (UBh)  +  k  <u2Bh)  +  <Vbb)k-m  -  ("b"bb)k-'i 


+  <PBh>  -  \  h  (UBh> +  -  (Tzb>k-.s  ■ 0 


(1) 


with 


of  the  surface  =  C* (pQ  /P)wa  cos  T 


of  an  interlayer  =  Az  9U  /  3z 


A 


X  of  the  bottom  =  g  U  U|/c 


/c 


•• 
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Cl 


continuity 


(w  t>)  , 

b  k--$ 


a  (Zb)  _  9_ 
9t  9  x 


(w,  b),  .  +  l-  ( LJBh ) 
b  k+S  9  x 


(UB)  dz  +  -  0 


_  2<2 
9  x 


< i nterna  1 ) 


vertical  (z-direction)  momentum 


3P 
9  z 


=  Pg 


heat  balance 


H  <BhT)  +  3J  (UBhT>  +  <”bbT)k-Hs  '  (“bbT>M* 


'<<_  (  SBh'r 
9x  I  x  9x 


salinity  balance 


9  BT 
9  z 


k+k 


9BT 


z  9z  J k-k 


H  Bh 
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<Bhs>  +  51  <UBhs>  +  <"bbs>k+ii  -  ("bbs,k-is 


{  9Bhs] 

f  9BS 

9BS 
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D  \ 

{  Z  dz  J 

H  Bh 
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constituent  balance 


Tt  +  (UBhc)  +  (wbbC)k+^  “  (wbbC)k-i-i 


i  9Rhc) 
Inx  Tx  J 


9BC 


k+k 


9BC 

z  ^ 


D  37" 


Bh 


J  k-k 


(2a) 


(over  total  depth)  (2b) 


(3) 


(4) 


(3) 


(6) 


* 


state 


1000  P 

P  =  _ o  . 

A  +  0. 6980  P0 


where 

A  =  1779.5  +  11.25T  -  0.0745T2  -  (3.80  +  0.01T)S 
P0  =  5890  +  38T  -  0.375T2  +  3S 

where 

A  x-direction  momentum  dispersion  coefficent  (m2/s) 

x 

A  z-direction  momentum  dispersion  coefficent  (in  /s) 

z 

b  estuary  width  (m) 

B  laterally  averaged  estuary  width  integrated  over  h  (m) 

% 

c  Chezy  resistance  coefficient,  m  /s 

C  laterally  averaged  constituent  concentration  integrated 

over  h  (mg*l-1) 

C*  resistance  coefficient 

x-direction  temperature  and  constituent  dispersion 
coefficient  (m2/s) 

0  z-direction  temperature  and  constituent  dispersion 

z 

coefficient  (m2/s) 

g  acceleration  due  to  gravity  (m/s  ) 

h  horizontal  layer  thickness  (m) 

H  total  depth  (m) 

H  source  strength  for  heat  balance  (C-m3*s~!),  salinity 

n 

(PPT *m3 • s-1 )  ,  or  consisttient  balance  (mg  •  1  ~ 1  •  m 3  •  s~ 1  ) 
k  integer  layer  number,  positive  downward 

P  pressure  (N/m2) 

Q  tributary  inflow  and  withdrawal  rates  (m'/s) 

salinity  (ppt) 


S 


T 

U 


u 


b 


V 


W 

a 


x  and 


laterally  averaged  temperature  integrated  over  h  (C) 
x-direction,  laterally  averaged  velocity  integra. _J  over 
h  (m/s) 

x-direction,  laterally  averaged  velocity  (m/s) 
cell  volume  (B'h’Ax)  (m3) 
wind  speed  (m/s) 

z-direction,  laterally  averaged  velocity  (m/s) 

Cartesian  coordinates:  x  is  along  the  estuary  centerline 
at  the  water  surface,  positive  downestuary,  and  z  is  posi¬ 
tive  downward  from  the  x-axis  (m) 
surface  elevation  (m) 
density  (kg/m3) 
air  density  (kg/m3) 
z-directicn  shear  stress  (m2/s2) 
wind  direction  (rad) 


8.  The  finite  difference  operation  is  applied  to  Equations  (1) 
through  (6)  and  introduces  two  additional  variables,  Ax  ,  the  x-direction 
spatial  step  (m)  and  At  ,  the  computation  time  step  (s) .  This  operation 
produces  finite  difference  equation  analogues  of  Equations  (1)  though  (6). 
The  solution  technique  is  to  substitute  the  forward  time,  longitudinal 
momentum  term  from  Equation  (1)  into  the  vertically  integrated  continuity 
expression,  Equation  (2b),  to  give  the  free  surface  frictionally  and  in- 
ertiallv  damped  longwave  equation.  The  latter  is  solved  implicitly  for 
the  surface  elevation  Z  ,  which  eliminates  the  longwave  speed  time  step 
restriction  (Ax/At  >  /gH)  .  However,  the  computation  time  step  At  is 
still  limited  by  the  fundamental  condition  that 


At  <  V/Q* 


(8) 


for  each  cell,  where  Q*  is  the  flow  into  or  out  of  the  cell  and  V  is 


9.  The  computation  proceeds  as  follows:  1)  new  water  surface  eleva¬ 
tions  are  computed  using  the  implicit  equation  combination  (Equations  1 
and  2b);  2)  new  longitudinal  velocities  on  the  grid  are  computed  from  the 
explicit  Equation  (1);  3)  Equation  (2a)  is  used  to  compute  vertical  veloc¬ 
ities  and  Equation  (2b)  is  used  to  check  water  balances  prior  to  the  im¬ 
plicit  solution  of  Equations  (4),  (5),  and  (6)  for  temperature  salinity 
and  constituent  concentrations;  4)  densities  are  computed  from 
Equation  (7)  (Leendertse  1973).  Results  are  printed,  and  a  new  time  step 
computation  begins  once  again  with  the  solution  for  water  surface 
elevation. 

10.  The  FORTRAN  names  of  the  variables  introduced  in  Equations  (1) 


through 

(7) 

are  given  below* : 

variable 

FORTRAN  name 

variable 

FORTRAN  name 

A 

AX 

S 

S1.S2 

X 

A 

AZ 

2 

b 

B 

T 

Tl  ,T2 

B 

B 

U 

U 

c 

CHZY 

Ub 

u 

C 

C1,C2 

V 

VOL 

C* 

cz 

w 

a 

WA 

D 

X 

DX 

Wb 

W 

D 

DZ 

Z 

Z1  ,Z2 

2 

g 

G 

At 

DLT 

h 

HIN.Hl ,H2 

Ax 

DLX 

tl 

HIN,H1,H2 

A 

LA 

H 

HN 

P 

RHO 

n 

k 

K 

Pa 

RHOA 

P 

P 

T 

z 

ST 

Po 

Q 

PO 

QTRIB,QWD 

0 

PHI 

*  A  glossary 
sented  in 

of  the  FORTRAN  variables 
Appendix  E. 

discussed  in  this 

report  is  pre- 

11.  The  location  of  the  variable  on  the  finite  difference  grid  is  im¬ 
portant  in  understanidng  the  FORTRAN  coding  of  CE-QUAL-ELV2  and  the  out¬ 
put.  Figure  1  shows  the  location  of  velocities  and  dispersion  coefficients 
at  cell  boundaries.  This  convention  permits  boundary  values  of  these  vari¬ 
ables  to  be  exactly  zero.  Other  variables  are  located  at  the  cell  center 
and  represent  averages  for  the  entire  cell. 

12.  Throughout  the  computation  new  values  of  variables  replace  old 
values.  The  solution  technique  requires  that  values  of  Z  ,  H  ,  T  , 

S  ,  and  C  be  retained  for  two  time  steps.  At  the  end  of  a  time  step, 
the  contents  of  these  variable  arrays  are  exchanged  so  that  the  older 
values  are  placed  in  arrays  with  a  "1"  suffix  and  the  newer  values  are 
placed  in  arrays  with  a  "2"  suffix. 

Estuarine  Boundary  Conditions 

13.  The  upstream  boundary  condition  for  an  estuary  is  the  same  as  for 
the  reservoir  case  and  consists  of  setting  left-hand  horizontal  velocities 
based  on  known  inflow  rate  and  geometry.  The  estuarine  or  open  boundary 
condition  is  more  complex.  Generally,  the  right-hand  boundary  water  sur¬ 
face  elevation  is  known  and  the  right-hand  horizontal  velocities  and  net 
flow  must  be  computed.  The  individual  components  of  the  x-momentum 
equation  are  computed  prior  to  the  solution  of  the  implicit  combination 
continuity  and  x-momentum  Equations  (1)  and  (2b).  This  solution  uses  the 
known  boundary  value  of  Z  ,  the  individual  components  of  the  x-momentum 
terms,  and  the  previous  values  of  Z  to  determine  new  values  of  Z  . 
Following  this  solution,  the  horizontal  velocities  U  are  explicitly 
computed . 

14.  Because  of  limitations  in  the  amount  of  information  available  be¬ 
yond  the  right-hand  boundary  (only  Z  and  the  temperature,  salinity,  and 
constituent  profiles,  not  velocities,  are  assumed),  several  of  the  compo¬ 
nents  in  Equation  (L)  cannot  be  computed  at  the  right-hand  boundary.  This 
equation  is  reduced  at  the  boundary  to: 
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rigure  1.  Location  of  major  variables  on  finite  difference  grid 


~  (UBh)  +  ^  (PBh)  +  (x  b),  ,,  -  (T  b),  , 
3 1  p  3x  z  k+'-s  z  k-!i 


Only  these  terms  at  the  right-hand  boundary  are  used  in  the  solution  for 
Z  and  U  .  The  net  inflow  or  outflow  at  this  boundary  is  then  computed 
by  integrating  the  U's  over  the  depth. 

15.  Temperature,  salinity,  and  constituent  profiles  are  also  needed 
at  the  boundary  in  order  to  solve  Equations  (4),  (5),  and  (6).  The 
right-hand  boundary  velocities  computed  from  Equation  (9)  are  used  to 
transport  these  constituents  into  or  out  of  the  estuary.  The  FORTRAN 
names  for  these  boundary  profiles  are  TR  ,  SR  ,  and  CR  ,  which  are  a 
function  of  depth.  The  FORTRAN  name  for  the  boundary  surface  elevation 


is  ZR 


Validation  Tests 


16.  The  codes  from  which  CE-QUAL-ELV2  was  derived  have  been  validated 
and  verified  (Edinger  and  Buchak  1979,  1980a,  and  1981).  The  validation 
tests  were  based  on  simulations  of  cases  for  which  consistency  and  accu¬ 
racy  of  water,  heat,  salinity,  and  constituent  balances  could  be  determined. 
Verification  tests  compared  field  observations  with  simulation  results. 

17.  Validation  tests  made  specifically  for  CE-QUAL-ELV2  focused  on  the 
newly-added  Water  Quality  Transport  Module.  These  tests  separately  exer¬ 
cised  features  of  the  program  that  handled  the  boundary  conditions  of  inflow, 
outflow,  and  exchange  at  the  open  water  boundary.  Each  test  was  deemed 
successfully  completed  if  perfect  temperature,  salinity,  and  constituent 
balances  were  obtained.  Evaluation  of  the  results  was  simplified  because 
each  of  the  cases  was  run  using  uniform  and  constant  temperatures,  salin¬ 
ities,  and  constituents,  with  no  sources  or  sinks  other  than  at  inflows, 
outflows,  or  the  open  boundary.  The  validation  tests  were  as  follows: 

Case  Descr  i  P  t_ion 


No  motion.  Constant  and  uniform  temperature,  salinity  and 
constituent;  no  inflow  or  withdrawals;  no  head,  tempera¬ 
ture,  or  salinity  gradient  at  open  boundary. 
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Case 


Description 

2  Inflow,  tributary,  and  withdrawal.  Same  as  1,  with  iso¬ 
thermal  and  isohaline  inflow,  tributary,  and  withdrawal. 

3  Surge.  Same  as  1,  with  open  boundary  surface  elevation 
higher  than  initial,  internal  elevations. 

4  Convection.  Same  as  1,  with  boundary  temperature  less 
than  initial,  internal  temperature  and  boundary  salinity 
greater  than  initial,  internal  salinity. 

5  Tide.  Same  as  1,  with  time-varying  open  water  boundary 
elevation. 


All  these  cases  were  successfully  run. 

18.  Figure  2  shows  the  displacement  vectors  derived  from  the  computed 
velocity  field  for  day  three  of  validation  case  4. 
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Figure  2.  Displacement  vectors  for  convec¬ 
tion  validation  case 


The  displacements  were  obtained  by  multiplying  the  velocities  by  the 
scaling  parameter  (in  this  case,  0.5  day)  to  obtain  distances  which  were 
then  plotted  on  the  finite  difference  grid.  The  figure  shows  movement  up 


the  estuary  In  the  bottom  and  movement  down  the  estuary  in  the  surface 
layers  in  response  to  the  density  gradient  at  the  open  water  boundary. 
Denser  boundary  water  (i.e.,  colder  and  more  saline)  entered  the  estuary 
in  the  lower  layers,  and  less  dense  estuary  water  left  in  the  surface 
layers.  This  process  continued  as  the  density  gradient  moved  up  the 
estuary  when  colder,  saltier  water  was  advected  in.  By  three  days  into 
the  simulation,  the  salinity  had  reached  segment  17  (km  32). 

19.  Validation  case  1  showed  that  extremely  small  velocities 
(<  l  x  10-8  m/s)  were  generated  at  a  large  number  of  time  steps  for  a 
case  in  which  zero  motion  was  expected.  These  spurious  velocities  were 
probably  due  to  the  accuracy  limits  of  the  computer  and  are  dependent 
on  the  size  of  the  grid  considered.  These  velocities  occurred  because 
of  a  surface  elevation  gradient  that  was  computed  as  the  result  of  many 
secondary  computations  according  to  the  recursive  algorithm  TRIDAG. 


PART  III:  LATERALLY  HOMOGENEOUS  ESTUARIES 


20.  Estuaries  have  been  classified  by  Pritchard  (1956)  on  the  basis 
of  their  dominant  salinity  gradients.  An  estuary  which  has  salinity 
variations  along  its  length,  over  its  depth,  and  across  its  width  is, 
of  course,  three-dimensional.  A  broad,  shallow  estuary  which  has  pre¬ 
dominately  longitudinal  and  lateral  salinity  gradients  but  no  vertical 
salinity  gradient  is  considered  a  vertically  homogeneous  estuary.  A 
long,  narrow,  and  relatively  deep  estuary  with  little  lateral  variation 
in  salinity  but  with  longitudinal  and  vertical  gradients  is  considered 
a  laterally  homogeneous  estuary  and  is  the  type  of  estuary  to  which  the 
CE-QUAL-ELV2  relationships  apply.  An  estuary  with  little  lateral  and  ver¬ 
tical  variation  in  salinity  but  with  a  longitudinal  gradient  is  considered 
to  be  a  sectionally  homogeneous  estuary. 

21.  The  above  classifications  provide  a  basis  for  simplifying  the 
hydrodynamic  and  transport  relationships  in  an  estuary  to  match  its 
spatial  variability.  One-dimensional  sectionally  homogeneous  hydro¬ 
dynamics  is  simpler  than  two-dimensional  vertically  homogeneous  hydro¬ 
dynamics,  which  in  turn  is  simpler  than  two-dimensional  laterally 
homogeneous  or  three-dimensional  hydrodynamics.  Complexity  of  the  re¬ 
lationships  is  determined  by  the  number  of  variables  being  computed, 

the  number  of  spatial  dimensions,  the  number  of  boundary  conditions,  and 
the  types  of  numerical  computations  that  must  be  employed.  The  numer¬ 
ical  hydrodynamics  of  all  three  classes  of  estuaries  has  been  developed 
by  Edinger  and  Buchak  (1980a). 

22.  Salinity  gradients  alone  may  not  dictate  the  proper  classifica¬ 
tion  of  an  estuary.  An  estuary  that  is  considered  sectionally  homoge¬ 
neous  based  on  salinity  can  exhibit  complex  time-varying  vertical 
velocity  profiles,  with  surface  velocities  reversing  before  bottom 
velocities  on  each  oscillatory  reversal  of  the  tide.  Use  of  sectionally 
homogeneous  relationships  in  this  case  results  in  mixing  because  the 
complex  velocity  shear  is  absorbed  into  the  evaluation  of  one-dimensional 


dispersion  coefficients.  Separate  consideration  of  mining  in  the  vertical 
and  longitudinal  axes  would  permit  more  realistic  evaluation  of  mixing 
coefficients . 

23.  Problems  of  oversimplifying  an  estuary  to  the  sectionally  homo¬ 
geneous  case  can  be  overcome  by  use  of  the  laterally  homogeneous  or  CE- 
QUAL-ELV2  relationships.  The  amount  of  boundary  data  (including  fresh 
water  inflows,  tide  heights,  and  tidal  salinity  at  the  downestuary  bound¬ 
er  > )  is  similar  for  the  two  cases.  The  CE-QUAL-ELV2  geometry  provides  a 
more  realistic  representation  of  cross  sections,  including  overbank  areas, 
and  a  more  accurate  location  of  intakes  and  discharges.  The  CE-QUAL-ELV2 
computations  cost  more  than  sectionally  homogeneous  relationships,  but  they 
overcome  the  problem  of  proper  representation  of  mixing  due  to  complex  ver¬ 
tical  velocity  structure,  and  they  provide  detail  about  constituents  that 
may  mix  slowly  along  the  vertical  axis. 

Dens  L_ty  Dynamics 

.’4 .  Many  estuaries  exhibit  a  complex  longitudinal  and  vertical  salin¬ 
ity  structure  due  to  tin;  intrusion  of  higher  density  salt  water  up  the 
bottom  and  under  the  lighter  freshwater  inflow  moving  downestuary.  There 
is  mixing  between  the  two  layers  due  to  tidal  action,  velocity  shear,  and 
wind  shear,  causing  the  salinity  structure  to  vary  within  tidal  cycles 
and  over  longer  periods  of  time.  The  vertical  velocity  structure  also 
varies  aid,  when  averaged  over  many  tidal  cycles,  shows  a  characteristic 
upHsttiary  flow  along  the  bottom  and  a  downestuary  flow  on  the  surface. 

2i.  The  characteristic  dynamics  of  a  laterally  homogeneous  estuary 
are  determined  hv  the  variation  of  the  horizontal  pressure  gradient  with 
depth.  1  he  influence  of  the  horizontal  density  gradient  on  the  hori- 
/. on t  a  1  pressure  gradient  can  be  demonstrated  algebracially  beginning 


with  tic  i ■  orient  it  inn  o 


pressure  P  from  a  free  surface  elevation 


to  ,a  depth  7  .  ibe  vertical  dis trihut  ion  of  pressure  is: 


where  p  is  the  density.  Since  p  and  p  vary  horizontally,  the 
horizontal  pressure  gradient  is: 


■z 
9P 

*  ■  -i 

The  horizontal  pressure  gradient  depends  on  the  surface  slope  9p/9x 
and  varies  vertically  with  the  integral  of  the  horizontal  density 
gradient.  For  purposes  of  demonstration,  assume  that  9p/9x  is  constant 
with  depth  as  would  be  the  case  for  a  sectionally  homogeneous  estuary. 

The  horizontal  pressure  gradient  then  beomes: 
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The  vertical  variation  of  the  horizontal  pressure  gradient  is  clear  from 
the  above  equation.  At  the  surface  z  =  p  ,  the  horizontal  pressure 
gradient  follows  the  surface  slope  as  9P/9x  =  gp(9p/9x)  .  At  an  inter¬ 
mediate  depth  of  z  =  p  +  (p3p/9x)y(9p/3x)  ,  the  horizontal  pressure 
gradient  is  zero  and  there  is  a  level  of  no  motion.  Below  this  depth, 
the  horizontal  pressure  gradient  is  opposite  in  direction  to  the  surface 
slope  due  to  the  vertical  integral  of  the  horizontal  density  gradient. 
The  level  of  no  motion  varies  with  time;  p  and  9p/9x  change  with 
propagation  of  the  surface  tidal  wave. 

26.  Even  a  sectionally  homogeneous  estuary  has  both  a  horizontal 
density  gradient  due  to  the  longitudinal  salinity  gradient  and  a  hori¬ 
zontal  pressure  gradient  that  varies  with  depth.  Velocities  in  a  sec¬ 
tionally  homogeneous  estuary  thus  vary  with  depth.  Over  a  tidal  cycle, 
surface  velocities  can  reverse  direction  before  bottom  velocities;  the 
net  result  is  a  two-layered  flow.  An  estuary  may  be  sectionally 


homogeneous  based  on  salinity  observations,  but  it  may  not  be  seetionally 
homogeneous  for  velocities  and  the  important  transport  processes.  A 
sectional  I v  homogeneous  estuary  based  on  salinity  would  be  well  repre¬ 
sented  by  laterally  averaged  estuarine  hydrodynamics.  Seetionally  homo¬ 
geneous  estuaries  are  a  special  case  of  laterally  averaged  dynamics. 

T 1  dal  Dynamics 

d'  t!,e  tide  is  one  of  the  main  driving  forces  of  an  estuary.  It 
the  surface  slope  of  the  estuary  to  which  the  internal  motions 
res  nond .  I'ldai  ierco  is  exerted  ori  the  estuary  in  the  rise  and  fall  of 
the  v.ito.r  surface  at  the  estuary  mouth  and  in  the  gravitational  pull  of 
the  moon  ana  sun  along  the  length  of  the  estuary.  The  latter  is  usually 
neglected  in  a  tidal  estuary  computation  because  the  estuary  is  domin¬ 
ated  by  the  propagation  of  the  tidal  wave  occurring  at  the  mouth. 

,.’b.  The  tidal  wave  at  the  mouth  propagates  up  the  esutary  at  the 
crv’itv  wave  speed  C  -  /gH  ,  which  varies  with  location  depending  on 
average  cross-sectional  depth.  For  an  estuary  with  an  average  depth  of 
}  m.  i a  .speed  of  7  m/s .  The  wave  length  L  is  dictated  by  the 
wave  period  "i  since  for  a  wave  C  =  L/T  .  For  a  tidal  period  of 
12.45  hours,  a  typical  wave  length  in  an  estuary  would  be  220  km.  For 
an  estuary  with  a  length  of  1 00  km  a  high  tide  may  exist  at  the  mouth 
while  a  low  1 1  tie  may  exist  at  the  head  of  tide.  A  typical  estuary  de- 

.  reuses  in  cross  sect  ion  and  depth  upestuary,  which  causes  a  tidal  wave 

t  s 1 . w  down  and  increase  in  amplitude  as  it  is  "f unneled"  through  the 

de>'  i a  -is  1  a-,  i.  ros  •;  sect  ion  . 

24.  Once  geometry  effects  are  accounted  for,  the  time  lag  of  a  tidal 

■  i'.c  peak  and  t. ho  increase  in  amplitude  of  the  tide  wave  upestuary  are 

r. :  1  led  1"'  b  a.  ton;  friction.  Accurate  tide  height  measurements  along 
i  i  )  f  •:  to  •■'■tu  i;---  give  the  tidal ly  averaged  surface  slope,  the 
:  i  iiou  it.i:T'e>mi  Lne  wave  time  Lai;.  All  three  are  used  to  adjust 

■  i  •.lome'-;.  rid  t  r  ici  u>n  since  they  dictate  the  characteristics  of 
t  i ‘e  1  i  d>.  w.'i iv. 

Wind  Dynamics 

i  •  iv!  1 . .  :  -in i  ir  an  osi-iary  depending  on  the  width,  fetch 

o  ■  I  ■;  > ;  i  :  -v  i  >>.•  1'stinrv  relative  to  prevailing  winds. 


Wind  effects  are  exerted  on  an  estuary  (a)  as  a  component  of  the  water 
surface  elevation  record  at  the  mouth  of  the  estuary  and  (b)  by  direct 
wind  shear  on  the  estuarine  water  surface.  Wind  effects  as  exhibited  by 
fluctuations  in  mean  water  levels  have  periods  ranging  from  hours  to  days. 

31.  When  determining  the  hydrodynamics  of  an  estuary,  the  wind  compo¬ 
nent  of  the  water  surface  elevation  recorded  at  the  mouth  and  the  surface 
shear  within  the  estuary  must  be  considered  together.  To  impose  wind  sur 
face  stress  within  the  estuary  and  tidal  height  variation  at  the  mouth  with¬ 
out  the  wind  component  would  produce  a  model  with  an  artificial  estuarine 
circulation  lacking  a  mass  transport  component  at  the  mouth.  Similarly,  to 
impose  an  observed  variation  in  water  surface  elevation  at  the  mouth  that 
does  contain  a  wind  component  without  considering  a  complimentary  wind 
surface  shear  would  produce  a  model  (a)  that  propagates  long-period  in¬ 
ternal  oscillations  not  compensated  for  by  surface  shear  and  (b)  with  an 
artifical  circulation. 

32.  For  analysis,  the  wind  component  is  often  filtered  out  of  the 
tide  record  at  the  estuary  mouth,  thus  eliminating  the  need  for  internal 
parameterization  of  wind  effects.  When  comparing  model  results  to 
periods  of  real  data,  however,  it  should  be  determined  if  the  data  were 
filtered.  Otherwise,  internal  parameters  may  be  erroneously  adjusted  to 
compensate  for  a  missing  driving  force. 

33.  Wind  plays  an  important  role  in  determining  the  motion  or  in¬ 
tensity  of  turbulence  at  spatial,  scales  shorter  than  the  Ax  or  Az 
scales  of  the  numerically  computed  mean  motions. 


Turbulent  Transport 

34.  The  source  of  turbulent  transport  terms  is  the  temporal  averaging 
of  the  hydrodynamic  and  transport  relationships;  temporal  averaging  gives 
time-smoothed  results  from  fluctuating  turbulent  flows.  The  turbulent 
transport  terms  are  derived  from  an  instantaneous  value  of  velocity  u 
represented  by  a  short-time  averaged  mean  value  U  and  the  fluctuation 
about  the  mean  u'  as: 


u  =  U  +  u' 


(13 


and  similarly  for  the  vertical  velocity  and  salinity: 

w  =  W  +  w'  (14 

s  =  S  +  s'  (15 

Averaging  the  horizontal  and  vertical  advection  of  momentum  terms 

~uiu/3x  and  9uw/9z  and  the  horizontal  and  vertical  advection  of  salinity 

9us/3x  and  9ws/9z  gives: 


uu  =  U  U  +  <u’  u’> 
uw  =  U  W  +  <u'  w’> 
ws=WS+<u'  s’> 


where  the  overbar  ( — )  signifies  a  time  average  of  the  product  terms.  The 
momentum  and  constituent  transport  balances  are  for  the  mean  values  U  , 

W  ,  and  S  ,  and  their  averagings  result  in  the  average  cross  product 
fluctuation  terms  <u'  u'>  ,  <u'  w'>  ,  <u'  s’>  ,  and  <w'  s'>  which  in 
turn  represent  the  horizontal  and  vertical  turbulent  fluxes  of  momentum 
.ind  salt. 

55.  A  description  of  the  turbulence  processes  in  a  numerical  model  is 
complicated  by  the  relationship  between  the  model  integration  time  step 
t  and  the  time  period  of  averaging.  Presumably,  the  mean  values  U  , 

W  .  and  S  are  an  average  over  the  integration  time  step,  and  the  fluc¬ 
tuations  about  the  mean  u'  ,  w'  ,  and  s'  are  measured  many  times 
within  this  time  step.  A  short  integration  time  step  of  a  few  seconds 
to  one  or  two  minutes  approaches  the  frequency  of  the  fluctuations,  yet 
tin.'  model  is  not  actually  computing  turbulent  fluctuating  velocities. 
Rather,  the  model  is  iterating  between  the  time  limits  of  the  boundary 
a  it  a  to  the  next  solution  point.  One  advantage  of  an  implicit  solution 
with  long  time  steps  is  1 1  in t  the  computed  U  ,  W  ,  and  S  at  the  end  of 
oil'll  time  step  represents  a  mean  over  a  period  for  which  the  average  of 
tin  fluctuations  apply. 
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36.  A  second  difficulty  in  describing  the  turbulent  transport  processes 
for  use  in  a  numerical  model  is  relating  them  to  the  scale  of  the  grid. 
Although  the  turbulent  transport  terms  are  properly  multiplied  by  the  width 
of  lateral  averaging,  the  mean  velocity  and  salinity  in  the  gradient  terms 
3U/3z  and  3S/3z  apply  across  a  segment  area  BAx  where  B  is  of  the 
order  of  hundreds  of  meters  and  Ax  is  of  kilometres.  The  width  average 
is  implied  when  computing  U  ,  W  ,  and  S  ,  but  the  length  average  is 
not.  Thus,  there  is  a  component  of  spatial  variation  contained  in 

<u'  w'>  and  <w'  s'>  that  varies  over  the  computational  grid  length  Ax  . 

Turbulence  Models 

37.  The  problem  of  relating  the  turbulent  transport  given  by  <u'  u'> 
and  <u'  w’>  for  momentum  and  <u'  s*>  and  <w'  s’>  for  constituent 
transport  to  the  mean  flow  field  quantities  of  U  ,  W  ,  and  S  is  known 
as  the  closure  problem.  There  are  two  general  classes  of  turbulent  closure 
models  that  have  been  studied  in  recent  years.  One  class  of  models  relates 
the  turbulent  transport  quantities  to  gradients  of  the  mean  flow  field 
through  dispersion  coefficients;  the  other  evaluates  the  turbulent  transport 
quantities  from  higher  order  moments  or  averages  of  the  momentum  and  trans¬ 
port  balances.  The  dispersion  coefficient  formulation  has  a  hierarchy  of 
models  that  have  developed  over  the  years  and  that  have  an  empirical  basis. 
The  higher  order  moments  models  have  been  applied  successfully  only  to  rela¬ 
tively  simple  problems  and  usually  require  an  empirical  closure  relation. 

Dispersion  Relations 

38.  The  dispersion  description  of  the  turbulent  transport  is  based  on 
an  analogy  to  molecular  diffusion.  The  turbulent  transports  are  related 
to  the  mean  flow  field  as: 
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The  turbulent  dispersion  coefficients  for  momentum  A  and  A  are  sume- 

x  z 

times  called  eddy  viscosities;  those  for  constituent  transport  W  and  U 
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are  sometimes  called  eddy  dif fusivities .  They  are  simply  proportionality 

coefficients  between  the  turbulent  transport  quantities  and  the  mean 

flow-field  quantities  and  are  not  implied  to  be  constant  in  time  or  space. 

39.  Rodi  (1980)  indicates  that  the  dispersion  coefficient  description 

of  the  turbulent  transport  processes  does  not  constitute  a  turbulence 

model  but  only  provides  the  framework  for  constructing  such  a  model,  and 

the  problem  shifts  to  evaluating  A  ,  A  ,  D  ,  and  D  over  space 

x  z  X  z 

and  time.  The  dispersion  coefficient  description  of  the  turbulent  trans¬ 
port  process  has  been  incorporated  into  the  LARM2  and  CE-QUAL-ELV2  hydro¬ 
dynamics  to  be  a  function  of  time  and  space  for  independent  evaluation  by 
any  formulation. 

AO.  There  is  a  hierarchy  of  models  for  evaluating  the  turbulent  dis¬ 
persion  coefficients  from  characteristics  of  the  mean  flow  field.  They 
are:  (a)  the  zero-equation  models  that  use  scaling  velocities  and 

lengths,  (b)  the  one-equation  models  that  utilize  the  transport  of  turbu¬ 
lent  kinetic  energy  in  the  flow  field  and  scaling  lengths  and,  (c)  the 
two-equation  models  that  utilize  the  transport  of  turbulent  kinetic  energy 
and  its  dissipation  in  the  flow  field.  The  latter  two  models  have  only 
recently  been  tested  for  application  in  detailed  estuarine  hydrodynamic 
models . 

Zero-Equation  Models 

41.  The  zero-equation  models  assume  that  the  turbulent  dispersion 
coefficients  are  proportional  to  a  scaling  velocity  V  and  scaling  length 

L  a  s : 

A  «  VL  (23 


which  shifts  the  burden  to  evaluating  V  and  L  .  The  scaling  length 
1,  characterizes  the  large-scale  turbulent  motion  related  to  the  mean 
flow  field.  In  finite  numerical  computations  its  evaluation  is  compli¬ 
cated  by  the  fact  that  the  mean  flow  field  is  computed  over  a  finite  Ax 


42.  The  PrandLl  mixing  length  hypothesis  relates  the  staling  velocity 
to  the  mean  flow  field  as: 


V  =  lm  |  3U/3z  |  (24) 

where  lm  is  a  mixing  length  related  to  boundary  layer  thickness. 

Pritchard  (1960)  has  used  a  mixing  length  hypothesis  to  evaluate  vertical 
momentum  dispersion  coefficients  frem  field  evaluations  in  the  James  River. 

43.  An  extension  of  the  Prandtl  mixing  length  hypothesis  is  von  Karman's 
evaluation  of  the  mixing  length  as: 

lm  =  k  |  (3U/3z)/(32U/3z2)  |  (25) 

which  applies  to  relatively  simple  flows.  Application  of  von  Karman's  re¬ 
lation  directly  to  estuaries  is  limited  by  the  fact  that  there  are  often 
zero  vertical  velocity  gradients  or  curvature  causing  lm  to  extend  from 
zero  to  infinity.  The  latter  requires  additional  truncation  hypotheses. 

44.  Zero-equation  formulations  of  the  dispersion  coefficients  are 
presently  utilized  in  CE-QUAL-ELV2  because  of  the  background  of  empirical 
information  available  for  their  evaluation  in  coastal  plain  estuaries. 

One-Equation  Models 

45.  The  turbulent  dispersion  coefficients  can  be  theoretically  related 
to  the  turbulent  kinetic  energy  of  the  mean  flow  field.  The  relationships 
take  the  form  of : 


A  =  Cu  K2  L  (26) 

where  Cu  is  an  empirical  constant,  K  is  the  turbulent  kinetic  energy, 
and  L  is,  as  previously,  a  scaling  length. 

46.  The  turbulent  kinetic  energy  is  considered  to  be  generated  and 


transported  by  the  mean  flow  and  dissipated  by  viscosity.  The  transport  of 
turbulent  kinetic  energy  in  laterally  averaged  flow  is,  as  for  any  scalar 
constituent : 


JBK  JLiBK  ,  3WBK 

- - +  — -  +  — — 
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SS(U,W,S) 


(27) 


where  SS(U,W,S)  are  the  sources  and  sinks  as  functions  of  the  mean  flow 
field.  The  transport  of  turbulent  kinetic  energy  can  be  incorporated  di¬ 
rectly  into  the  transport  module  of  the  LAEM  code. 

47.  Leendertse  and  Liu  (1977)  have  utilized  the  one-equation  formula¬ 
tion  of  turbulent  dispersion  in  a  three-dimensional  time-varying  finite 
difference  computation  with  some  success.  The  formulation  allows  evalua¬ 
tion  of  the  generation  and  decay  of  turbulent  kinetic  energy  by  wind 
shear  at  the  water  surface.  The  buoyancy  component  of  the  source-sink 
term  allows  evaluation  of  the  dampening  of  turbulent  kinetic  energy  under 
stratified  conditions  and  the  conversion  of  potential  to  kinetic  energy  by 
unstable  density  gradients. 

48.  The  one-equation  or  K  model  requires  an  evaluation  of  a  scale 

length  for  the  particular  flow  problem  being  evaluated.  The  dissipation 
portion  of  the  source-sink  term  is  usually  modeled  as  £  =  K  '  /L 
where  is  another  empirical  constant.  Often  in  simple  flows,  the 

length  scale  is  evaluated  from  the  von  Karman  formulation. 

f wo-Equation  Models 

49.  Some  of  the  empiricism  of  the  one-equation  or  K  model,  such  as 
scale  lengths,  can  be  removed  by  transporting  K  and  £  as  separate 
quantities,  each  with  a  separate  mean  flow  transport.  Both  quantities  can 
be  formulated  in  the  transport  module  of  CE-QUAL-ELV2 .  The  dispersion  co¬ 
efficient  then  becomes  a  function  of  both  K  and  z  as: 

A  =  Cu  K2/£  (28) 

and  the  length  scale  is  eliminated. 

50.  Rodi  (1980)  points  out  that  an  exact  dissipation  or  £  equation 
contains  complex  correlations  whose  behavior  is  little  known  and  for  which 
drastic  modeling  assumptions  must  be  made.  Carter  (1981)  has  reported 
that  numerical  experiments  by  Blumberg  in  atmospheric  stratified 
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flows  using  the  K-e  relationships  did  not  Improve  flow-field  predic¬ 
tions  much  beyond  those  given  by  the  zero-equation  models.  The  K-e 
relations  are  still  in  the  experimental  stage. 


Moment  Models 

51.  It  is  possible  to  formulate  relationships  for  the  turbulent  trans¬ 
port  terms  <u'  u’>  ,  <u'  w'>  ,  <u*  s’>  ,  and  <w'  s'>  directly  from 

the  momentum  and  transport  equations.  For  example,  multiplying  the  salin¬ 
ity  transport  by  u  =  U  +  u'  and  averaging  gives  a  relationship  of  the 
form: 

3<u'  s'>  ,  3U<u'  s’>  3<u'  u'  s'>  . 

- - -  + - - - - -  +  •••  (29 

3t  3x  3x 

which  introduces  higher  moment  terms  such  as  <u'  u'  s’>  ,  <u'  w'  s’>  » 

and  so  on.  Thus,  Introducing  a  relationship  for  direct  evaluation  of 
<u*  s'>  results  in  higher  moments  that  need  to  be  evaluated.  The  process 
can  be  continued  to  fourth-order  and  higher  moments.  However,  eventual 
solution  or  closure  requires  truncation  of  the  processes  and  empirical 
evaluation  of  the  highest  moment. 

52.  The  higher  moment  relationships  are  often  truncated  by  assuming 
that  the  turbulence  has  some  mathematical  property  such  as  isotropy.  The 
relationships  have  been  evaluated  for  simple  cases  such  as  flows  through 
grids  but  have  yet  to  be  applied  in  complex  time-varying  flows  such  as 
estuaries.  Part  of  the  reason  is  that  more  computations  are  required  for 
the  turbulent  transports  than  for  the  mean  flows. 

Dispersion  Relationships  in  CE-QUAL-ELV2 

53.  The  turbulent  transport  relationships  in  CE-QUAL-KLV2  follow  the 
zero-equation  models  and  use  dispersion  relationships.  The  dispersion 
coefficients  are  written  into  all  numerical  expressions  for  general  eval¬ 
uation  as  functions  of  time  and  space.  They  can  be  evaluated  from  scaling 
relationships,  as  shown  below,  or  they  can  eventually  be  evaluated  from  a 
K  or  K-f.  model.  the  scalar  relationships  rely  on  empirical  evaluation 
from  previous  estuary  studies  and,  where  possible,  field  verification. 


Vertical  Transport  Processes 

54.  The  average  cross-product  fluctuation  terms  which  moke  up  t lie 
vertical  turbulent  process  are  related  to  the  mean  flow  and  concentration 
fields  as: 


<u'  w'>  =  -A  3U/3z 
z 


<w'  s’> 


-D  3S/3z 


(30) 

(31) 
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The  problem,  therefore,  is  reduced  to  evaluating  and  Dz  in  terms  of 

the  computed  flow  and  density  fields. 

55.  In  an  estuary  that  can  stratify,  the  vertical  dispersion  coeffi¬ 
cients  are  affected  by  the  level  of  turbulence  and  the  degree  of 
stratification.  They  can  be  represented  as: 


t 

» 

t 


A  =  A  F (Ri)  (32) 

z  zo 

where  A  is  a  function  describing  the  level  of  turbulence  related  to 
zo 

velocity,  depth,  wind  stress,  and  Ri  is  the  Richardson  number  related  to 
stratification.  Almost  every  formulation  of  A^  that  is  available  can  be 
examined  in  the  above  form. 

56.  The  Richardson  number  is  an  important  concept  in  describing  verti¬ 
cal  mixing  in  a  stratified  waterbody.  It  is  defined  as: 
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which  is  the  ratio  of  the  buoyant  or  potential  energy  to  the  kinetic 
energy  being  dissipated  at  a  point  in  the  water  column.  This  Ri  is 
called  a  "local"  or  gradient  Richardson  number  as  opposed  to  a  "bulk" 
Richardson  number  applied  to  the  whole  water  column.  The  stronger  the 
stratification,  as  indicated  by  a  large  density  gradient  3p/3z,  the 
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more  a  given  amount  of  kinetic  energy  OU/9z)z  is  dissipated  by  the 

buoyancy  and  less  by  the  turbulence  generated.  An  increasing  Ri  , 

therefore,  means  a  lower  level  of  turbulence  and  a  lower  A  and  D 

7.  z 

57.  In  the  absence  of  stratification,  Ri  =  0  and  A  =  A  .  These 

v.  zo 

relationships  as  well  as  dimensional  arguments  led  Munk  and  Anderson  (1948) 
to  the  Richardson  number  function  of: 


F(Ri)  = 


(34) 


where  they  deduced  for  vertical  momentum  transport  a  =  10  and  b  =  -1/2 

and  for  vertical  salinity  transport  a  =  10/3  and  b  =  -3/2  .  Since  the 

exponent  b  is  negative,  a  higher  Ri  leads  to  a  lower  A  .  The  form 

z 

of  the  Richardson  number  function  given  above  appears  in  many  studies. 

It  tends  to  have  a  theoretical  basis,  as  shown  by  the  heuristic  deriva¬ 
tions  given  by  Officer  (1976). 

58.  Another  form  of  the  Richardson  number  function  used  by  Leendertse 
and  Liu  (1975)  is: 


F (Ri)  <=  e~r(R1)  (35) 

which  is  an  exponential  form  that  has  a  value  for  negative  Ri  and  can 
degenerate  to  zero.  A  negative  Ri  exists  where  an  inverse  density 
gradient  occurs  and  mixing  takes  place  by  turnover  or  by  vertical  con¬ 
vection  at  a  very  rapid  rate.  It  can  be  handled  in  numerical  computations 

by  using  an  algorithm  that  sets  a  large  A  where  the  density  gradient 

z 

is  negative.  A  large  Ri  can  lead  to  small  ,  the  lower  limit  of 

which  must  be  the  molecular  diffusivity  and  viscosity. 

59.  The  most  complete  evaluation  of  A  and  F(Ri)  for  an  estuary 

zo 

has  been  given  by  Pritchard  (1960) .  Based  on  data  from  the  James  River 
and  using  mixing  length  arguments,  he  deduced  that: 
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F (Ri)  =  (1  +  0.276  Ri) 


(36) 
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where  U  is  the  RMS  tidal  velocity  and  d  is  the  total  depth.  The 

Richardson  number  was  defined  as  (g/p  •  9p/3z)/(0. 70U/h)  where  U  is 

the  mean  velocity  over  the  water  column.  The  function  A  is  zero  at 

zo 

the  surface  and  bottom  of  the  water  column  and  has  a  maximum  at  middepth 
based  on  mixing  length  arguments.  This  shape  is  modified  by  F(Ri) 
which  reduces  the  middepth  value  because  of  stratification.  The  Ri  is 
determined  from  the  local  gradient  and  bulk  (water  column)  velocity  and 
depth.  It  is  clearly  intended  to  be  a  local  or  depth-dependent  evalua¬ 
tion  of  A^  by  choice  of  the  depth  function. 

60.  The  depth-dependent  forms  of  A were  tried  in  numerical  models 
by  Bowden  and  Hamilton  (1975)  and  by  Elliott  (1976a).  Both  reported 
numerical  instabilities  that  were  apparently  related  to  the  evaluation 
of  Az  using  a  Richardson  number.  Bowden  and  Hamilton  (1975)  resorted 
to  using  a  bulk  Richardson  number,  but  with  a  vertical-shaped  function. 
More  recently,  Bowden  (1977)  has  discussed  the  difficulties  of  utilizing 
the  local  Richardson  number  in  a  numerical  model  and  suggests  that  there 
is  a  theoretical  rationale  for  the  bulk  Ri  .  Examination  of  the  numer¬ 
ical  form  of  the  Elliott  and  the  Bowden  and  Hamilton  models  suggests  that 
the  instabilities  are  a  result  of  computational  procedures  rather  than 
physical  aspects  of  the  mixing  problem. 

61.  CE-QUAL-ELV2  presently  uses  a  local  Ri  determined  from  the 
local  velocity  and  density  gradient.  When  computations  are  made  for  an 
estuary  with  a  constant  Azq  ,  an  A^  is  produced  with  a  depth  distri¬ 
bution  as  given  by  Pritchard  (1960),  This  suggests  that  the  distribu¬ 
tion  of  A^  is  a  result  of  the  interaction  of  the  mean  velocity  and 
salinity  fields  with  the  vertical  dispersion. 

62.  Wind  is  important  in  vertical  mixing.  Wind  not  only  increases  the 
mean  velocity  gradient  and  shear,  but  must  also  increase  the  turbulent 
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fluctuations  about  the  mean.  It  is  the  fluctuating  component  over  a  given 
period  of  averaging  that  determines  the  vertical  mixing.  Pritchard  (1960) 
assumed  that  the  wind-induced  turbulence  is  proportional  to  the  orbital 
wave  velocity  resulting  from  the  wind  and  is  also  affected  by  the  vertical 
density  gradient.  The  wind  speed  function  deduced  by  Pritchard  (1960) 
from  the  James  River  data  is: 


A  =  9.57  x 10~3 

z(h  -  z) 

H)  e"2lIZ/'L 

zo 

h 

W  J 

(38) 


where  H  is  the  wind  wave  height,  T  is  the  wave  period,  and  L  is  the 
wave  length.  In  the  above  relationship,  the  Azq  decreases  exponentially 
with  depth  and  shows  a  parabolic  shape  with  depth.  The  exponential  decay 
occurs  because  the  wave  orbital  velocity  decays  with  depth,  as  derived 
from  elementary  wave  theory;  the  parabolic-shaped  terms  come  from  the 
mixing  length  theory  used  to  scale  to  the  total  water  column  depth. 

Ford  (1976)  has  examined  the  wind  effects  on  k^  primarily  for  lakes 
using  the  exponential  decay,  but  no  shape  function  with  depth.  The 
wind  wave  characteristics  H  ,  T  ,  and  L  are  functions  of  wind  speed, 
duration,  and  fetch.  The  functional  relationships  among  these  variables 
are  known  to  a  limited  extent  for  estuaries.  Accurate  representation  of 
vertical  mixing  in  estuaries  will  require  performing  wind  wave  analyses 
for  significant  wind  events.  The  scale  of  the  computational  model  grid 
relative  to  the  waterbody  has  a  significant  effect  on  the  magnitude  of 
the  dispersion  coefficients.  Increasing  layer  thickness  from  0.3  m  to 
1  or  2  m  reduces  the  resolution  of  the  computations;  the  reduced  resolution 
may  mean  larger  values  for  k^Q  .  In  principle,  when  the  grid  size  de¬ 
creases,  values  of  A  approach  molecular  scales.  However,  the  turbu- 

zo 

lence  in  the  waterbody  consists  of  random  motions  and  still  requires 
evaluation  of  the  turbulent  transport  terms  <u'  w'>  and  <w'  s'>  to 
close  the  equations  of  motion  and  transport.  In  a  unique  experiment  in 
model  scales,  CE-QUAL-ELV2  computations  were  performed  on  a  very  small  grid 


where  Ax  =  5  ft  and  H  =  3  in  a  hydraulic  flume  with  a  density  under¬ 
flow  at  laminar  conditions  (Edinger  and  Buchak  1979) .  The  numerical 
model  did  not  produce  accurate  results  until  the  dispersion  coefficients 
were  reduced  to  molecular  values,  indicating  that  there  is  a  strong 
interrelationship  between  grid  size  and  the  dispersion  coefficients  and 
turbulent  characteristics  of  the  waterbody. 

Longitudinal  Dispersion  Formulations 

63.  The  longitudinal  dispersion  of  momentum  and  constituent  concentra¬ 
tions  results  from  the  time  averaging  that  produces  the  advection  of 
momentum  9UU/9x  in  the  momentum  balance  and  the  constituent  9US/9x  in 
the  transport  balance.  The  average  of  the  product  fluctuations  about  the 
mean  values  are  for  momentum: 


<u '  u ' >  =  -A  9U/9x 
x 


and  for  constituent 


(39) 


<u '  s’>  =  -I)  9S/9x  (AO) 

x 

which  relates  the  turbulent  fluctuation  to  the  average  velocity  U  and 
constituent  S  .  As  discussed  previously,  there  is  a  relationship  between 
tne  time  over  which  the  fluctuation  products  are  averaged  and  the  time 
step  of  numerical  integration  that  produces  U  and  S  . 

64.  Most  formulations  of  longitudinal  dispersion  in  estuaries  have 
been  derived  in  relation  to  one-dimensional  sectionally  homogeneous  estuary 
models.  They  often  come  from  analysis  of  steady  channel  flows  and  are 
dimensionally  of  the  form: 


D  =  a  h*u* 
x 


(Al) 


where  h*  and  u*  are  a  characteristic  depth  and  velocity.  Harleman 
(1971)  shows  that  when  the  Taylor  formula  is  converted  from  pipe  flow  to 


a  tidal  case,  the  characteristic  velocity  translates  from  the  boundary 
shear  velocity  to  the  maximum  tidal  velocity.  Examples  of  the  above  rela¬ 
tionships  are  given  in  Fischer  et  al.  (1979). 

65.  Numerous  evaluations  of  have  been  made  for  specific  estuaries 

using  a  one-dimensional  model  of  the  salinity  distribution  (Officer  1977). 
For  the  Potomac  estuary,  values  of  10  to  20  m2/s  were  found  in  the  vicinity 
of  Washington,  D.C.,  increasing  to  60  m2/s  toward  the  mouth.  Using  a 
dissolved  oxygen  model  on  the  Delaware  River,  values  ranged  from  120  to 

210  m2/s  over  135  km.  The  Thames  exhibited  values  between  50  and  80  m2/s 
at  low  river  flows  and  up  to  340  m2/s  at  high  river  flows.  The  values, 
therefore,  can  range  from  5  to  500  m2/s  and  vary  with  river  flow  and 
salinity  stratification.  The  latter  can  complicate  results  since  the 
large  required  in  a  one-dimensional  model  in  the  stratification 

region  simply  forces  the  model  to  fit  a  multilayered  circulation,  and  a 
similarly  large  would  not  be  required  for  a  laterally  averaged  model 

in  the  same  region. 

66.  There  have  been  few  attempts  to  derive  relationships  for  the  use 

of  and  in  laterally  averaged  hydrodynamics  except  to  use  rela¬ 

tionships  as  if  they  applied  to  each  layer.  Experience  with  CE-QUAL-ELV2 
has  shown  that  changes  in  A^  and  of  two  orders  of  magnitude  have 

not  changed  results  significantly.  Similar  conclusions  were  reached  by 
Elliott  (1976a).  The  reason  for  this  lack  of  sensitivity  is  that  the 
contribution  of  vertical  exchange  to  longitudinal  mixing  which  is  absorbed 
in  in  one-dimensional  models  is  explicitly  included  in  the  laterally 

averaged  models.  Thus,  one  of  the  spatial  dimensions  over  which  the 
turbulent  transport  is  averaged  in  the  one-dimensional  model  is  relaxed 
in  the  two-dimensional  model  and  becomes  less  important. 

67.  There  still  may  be  a  scale  relationship  to  the  grid  Ax  for  A 

and  as  discussed  previously  for  the  flume  studies  at  laminar  flow 

conditions.  Results  suggest  that  in  each  computational  grid  it  is 
necessary  to  satisfy  the  condition  of: 
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Dmolecul nr  <  l)x  <  UAx 


(42) 


which  dimensionally  states  that  a  constituent  (including  momentum)  should 

not  be  diffused  out  of  a  cell  faster  than  it  is  advected  inward.  For  an 

estuary  like  the  Potomac  where,  computationally,  Ax  might  be  9300  m 

(5  nm)  and  the  average  net  non-tidal  velocity  is  of  the  order  of  0.02  m/s, 

maximum  D  would  be  186  m2/s. 
x 

Estuarine  Characteristics 

68.  A  laterally  homogeneous  estuary  is  usually  identified  by  its 
geography,  hydrography,  and  salinity  distributions.  A  laterally  homo¬ 
geneous  estuary  is  typically  defined  as  exhibiting  longitudinal  and 
vertical  gradients  and  lateral  homogeneity  of  salinity.  Velocity  pro¬ 
files  may  provide  a  more  accurate  basis  for  a  dimensional  description 
since  an  estuary  may  have  vertically  homogeneous  salinity  yet  exhibit 
velocity  reversals  with  depth.  Salinity  profiles  are  used  as  indicators 
because  of  the  availability  of  salinity  data.  A  laterally  homogeneous 
estuary  is  characteristically  a  coastal  plain  estuary  at  the  mouth  of  a 
drowned  river  valley.  The  estuary  is  much  longer  than  it  is  wide,  although 
there  may  be  a  bay  at  its  mouth.  Depths  may  be  shallow  and  cross-sectional 
areas  small  near  the  head  of  tide,  with  both  increasing  in  size  toward 

the  mouth. 

69.  Along  the  east  coast  of  the  United  States,  typical  laterally  homo¬ 
geneous  estuaries  are  the  Connecticut  River,  the  Hudson  River,  most  of  the 
tributaries  to  the  Chesapeake  Bay  including  the  Severn,  Patuxent,  Potomac, 
and  James  Rivers.  Long  Island  Sound  and  Chesapeake  Bay  can  be  approximated 
as  laterally  homogeneous  estuaries,  although  there  are  significant  lateral 
salinity  and  velocity  gradients  due  to  Coriolis  acceleration.  Many  of  the 
rivers  along  the  North  and  South  Carolina  and  Georgia  coasts  also  have 
sections  of  estuary  that  are  laterally  homogeneous  although  they  are  tribu¬ 
tary  to  shallow-water,  vertically  homogeneous  coastal  regions. 
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The  Patuxent  Estuary 

70.  CE-QUAL-ELV2  data  requirements  for  simulation  of  a  laterally 
homogeneous  estuary  are  best  illustrated  by  an  example.  The  Patuxent 
Estuary  was  chosen  because  it  shows  a  complete  transition  from  a  salinity- 
stratified  laterally  homogeneous  estuary  to  a  shallow,  sectionally 
homogeneous  estuary;  it  also  has  a  single  large  heat  source  which  allows 
the  demonstration  of  the  Water  Quality  Transport  Module  for  excess  heat. 

71.  The  Patuxent  Estuary  is  shown  in  Figure  3.  It  extends  70  km  from 
Solomons  Island  to  Hardesty,  Md.  Detailed  hydrographic  data  is  available 
for  the  Patuxent  in  numerical  form  (Cronin  and  Pritchard  1975)  .  Fresh¬ 
water  inflows  vary  widely  through  the  year,  with  high  flows  in  the  spring 
and  low  flows  in  the  fall.  Between  1970  and  197S  mean  monthly  flows  have 
varied  from  73  m3/s  to  2  m3/s.  Daily  freshwater  inflows  at  the  head  of 
the  estuary  and  at  major  tributaries  can  be  determined  from  the  United 
States  Geological  Survey  gage  records. 

72.  The  Patuxent  Estuary  is  tributary  to  the  Chesapeake  Bay.  The 
downestuarv  boundary  salinity  is  determined  by  conditions  in  the  bay. 

The  temporal  variation  and  stratification  of  salinity  at  the  mouth  of  the 
Patuxent  is  due  to  processes  within  the  Chesapeake  Bay,  including  the  in¬ 
fluence  of  the  Susquehanna  River. 

73.  The  salinity  distributions  within  the  Patuxent  River  estuary  are 

dictated  by  the  freshwater  runoff,  the  estuarine  geometry,  and  the  salin¬ 
ity  distribution  at  the  mouth.  Observed  longitudinal  and  vertical  dis¬ 
tributions  of  salinity  are  given  for  high  freshwater  inflow  period  and  a 
low  flow  period  in  Tables  1  and  2.  both  cases,  as  shown  by  the  verti¬ 

cal  salinity  profiles,  the  estuary  is  stratified  at  its  mouth  and  signif¬ 
icantly  through  its  lower  half.  In  the  headwater  reaches,  the  salinity 
is  significantly  lower  during  high-flow  periods  than  during  low-flow 
periods . 

lb.  The  spatial  and  temporal  vertical  distributions  of  salinity  at 
the  estuary  mouth  are  necessary  boundary  conditions  and  should  be  recorded 
at  daily  intervals,  if  possible.  Daily  surface  salinity  data  are  avail¬ 


able  for  the  Patuxent  at  Solomons  Island.  For  a  few  years  there  have 


Table  1.  Longitudinal  and  vertical  salinity  profiles  in  the  Patuxent  River  for  1?  May  1978,  PPT 


existed  monthly  vertical  salinity  profiles  that  must  be  interpolated  to 
daily  values.  Such  an  interpolated  boundary  record  presents  the  possi¬ 
bility  that  large,  short-term  fluctuations  in  salinity  due  to  storm 
events  will  be  missing  from  the  boundary  data  and  will  not  be  reflected 
in  the  computed  results. 

75.  Tide  height  data  at  the  estuary  mouth  are  a  necessary  boundary 
condition.  Tide  height  measurements  at  the  estuary  mouth  should  be  made 
accurately  and  taken  hourly,  if  possible.  Hourly  records  should  include 
the  effects  due  to  wind  setup  as  well  as  the  soiilunar  harmonics.  For 
the  Patuxent  at  Solomons  Island,  only  daily  high  and  low  tide  heights 
were  available,  and  many  data  were  missing  from  the  tide  record.  There¬ 
fore,  tide  heights  at  Solomons  Island  were  reconstructed  for  a  simulation 
period  using  the  following  data  sources:  observed  tide  heights  upestuary 
at  Benedict,  the  National  Oceanic  and  Atmospheric  Administration  (NOAA) 
predicted  tidal  amplitudes  for  the  estuary  at  Benedict  and  Solomons 
Island,  and  an  assumed  sinusoidal  tide  at  the  mouth.  Deviations  of  ob¬ 
served  and  predicted  tidal  amplitudes  at  Benedict  were  superimposed  on 
predicted  tidal  amplitudes  at  Solomons  Island  in  order  to  create  a  tidal 
boundary  record.  The  reconstructed  tide  record  at  Solomons  Island  does 
not  accurately  reflect  the  effects  due  to  wind  setup  either  along  the 
Chesapeake  Bay  or  on  the  estuary.  Tidal  amplitudes  are  among  the  most 
important  types  of  boundary  data  for  hydrodynamic  simulations,  yet  for 
most  estuarine  applications  tide  height  data  at  the  mouth  is  sparse  or 
nonexistent.  Tt  must  be  emphasized  that  the  accuracy  of  a  hydrodynamic 
model  is  no  greater  than  that  of  the  boundary  data  used. 

76.  A  major  heat  source  on  the  Patuxent  located  58  km  upestuary  from 
the  mouth  is  the  Chalk  Point  Power  Plant,  which  withdraws  high-salinity 
water  downestuary  and  then  discharges  it  upestuary  at  a  temperature  in¬ 
crease  of  7  C.  The  plant  pumps  at  a  rate  of  approximately  47  m3/s, 
representing  nnnv  times  the  fall  low  flow  and  a  large  fraction  of  the 
springtime  freshwater  flow.  The  more  saline  discharged  waters  are 
warmer  and  less  dense  than  tile  receiving  waters,  so  they  float  on  the 


surface  near  the  discharge.  As  they  mix  and  cool,  however,  their  density 
increases  relative  to  the  receiving  waters  and  they  sink.  These  types 
of  details  should  be  part  of  the  boundary  data  for  a  limited  region  of 
an  estuary. 

77.  There  have  been  a  number  of  intensive  surveys  of  the  Patuxent 
River,  including  the  monthly  salinity  and  water  quality  profiling  pre¬ 
sented  above,  that  provide  data  for  verification  of  simulations.  Velocity 
comparisons  present  one  of  the  most  stringent  verifications.  In  August 
1978,  extensive  water  current  surveys  were  made  in  a  region  that  extended 
between  35  km  and  45  km  from  the  estuary  mouth  (Academy  of  Natural  Sciences 
of  Philadelphia  and  J.  E.  Edinger  Associates,  Inc.  1981).  At  one  partic¬ 
ular  section  of  interest  at  Benedict,  four  surface  current  meters  and  four 
bottom  current  meters  placed  across  the  estuary  made  it  possible  to  deter¬ 
mine  laterally  averaged  velocities  at  two  depths  and  compare  them  to  com¬ 
puted  results. 

78.  Another  verification  compares  longitudinal  variations  in  tidal 
amplitude  and  mean  tide  height.  Except  for  the  Benedict  gage  these  are 
known  from  the  NOAA-predicted  tide  ranges  along  the  estuary,  and  the  com¬ 
parison  is  between  two  different  methods  of  computation.  Salinity  pro¬ 
files  within  the  period  can  also  be  compared;  either  long  time  intervals 
(30-60  days)  or  a  unique  storm  event  was  required  to  produce  a  uniquely 
different  salinity  profile.  The  Patuxent  simulations  were  initialized 
from  an  observed  salinity  distribution. 


PART  IV:  APPLICATION  PROCEDURE 
AND  EXAMPLE 


79.  Applying  CE-QUAL-ELV2  to  a  particular  case  requires  five  steps: 

(1)  assembling  the  necessary  data,  (2)  schematizing  the  estuary  geometry, 

(3)  modifying  a  subroutine  for  the  input  of  time-varying  boundary  condi¬ 
tion  data,  (4)  determining  site-dependent  hydraulic  constants  and  con¬ 
stituent  reaction  rates,  and  (5)  executing  the  code  and  evaluating  the 
results.  These  steps  are  described  below 

Required  Data 

80.  The  initial  task  in  applying  CE-QUAL-ELV2  to  a  specific  case 
(after  deciding  that  the  problem  requires  a  two-dimensional,  time-varying 
simulation  for  its  solution)  is  assembling  the  necessary  data.  Estuary 
geometry  is  of  immediate  usefulness  and  should  include  bathymetric 

cross  sections  (y-z  plane).  A  plan  (x-y  plane),  an  elevation  (x-z  plane), 
and  an  elevation-area-volume  table  are  useful,  but  not  necessary, 
adjunct  information.  These  data  are  used  to  construct  the  computational 
grid. 

81.  Initial  and  boundary  condition  data  are  required.  The  former 
are  measurements  of  water  surface  elevation  and  temperature  and  salinity 
fields  for  the  starting  day  of  the  simulation.  Boundary  condition  data 
in  their  most  general  form  consist  of  time-varying  inflow  rates,  tempera¬ 
tures,  and  salinities  for  the  upstream  inflow  and  each  significant  tribu¬ 
tary;  time-varying  release  and  withdrawal  rates;  and  time-varying  meteo¬ 
rological  data  for  wind  shear,  evaporation,  and  surface  heat  exchange 
calculations.  As  CE-QUAL-ELV2  is  now  coded,  surface  heat  exchange  is 
computed  using  equilibrium  temperature,  the  coefficient  of  surface  heat 
exchange,  and  solar  radiation.  These  parameters  are  computed  separately 
from  air  temperature;  dew  point  temperature;  wind  speed;  and  observed  solar 
radiation,  percent  of  possible  sunshine,  or  cloud  cover  data  (Edinger, 
Brady,  and  Geyer  1974).  In  their  simplest  form,  the  boundary  condition 
data  are  constants.  For  each  water  quality  constituent  to  be  computed, 
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initial  and  boundary  condition  data  and  reaction  rates  are  needed.  Hydrau¬ 
lic  properties — Chezy  and  dispersion  coefficients — and  solar  radiation  ab¬ 
sorption  and  attenuation  characteristics  are  also  required. 

Geometric  Schematization 

82.  The  geometric  data  are  organized  into  a  grid  like  the  one  shown  in 
Figure  4  for  the  Patuxent  Estuary.  The  basic  parameters  to  be  selected 

in  defining  the  grid  are  the  longitudinal  spacing  Ax  (FORTRAN  DLX)  in 
meters  and  the  vertical  spacing  h  (FORTRAN  HIN)  in  meters.  These  para¬ 
meters  are  constant  throughout  the  grid  and  define  two  of  the  three  di¬ 
mensions  of  each  cell.  The  third  dimension,  the  cell  widths  (FORTRAN  B) , 
are  most  easily  taken  from  cross-section  drawings  and  represent  an 
average  width  for  each  cell.  The  vertical  columns  defined  by  Ax  are 
called  segments,  and  the  horizontal  rows  defined  by  h  are  called  layers. 

83.  A  more  sophisticated  procedure  for  obtaining  cell  widths  involves 
the  use  of  the  Hydrologic  Engineering  Center  code  GEDA  and  the  prepro¬ 
cessor  GIN,  the  user  notes  for  which  are  included  as  Appendix  F  of  this 
guide.  GEDA  utilizes  cross-sectional  data  from  ranges  at  irregular  In¬ 
tervals  and  produces  estuary  widths  as  functions  of  elevation  at  regular 
intervals.  GEDA  also  computes  areas  and  volumes  as  functions  of  elevation. 
Cross-sectional  data  may  be  modified  using  GEDA  procedures  in  order  to  fit 
the  computed  elevation-area-volume  table  to  published  information.  GEDA 
allows  several  Ax's  and  h's  to  be  tried  with  very  little  additional 
effort.  GIN  uses  GEDA  output  as  input  and  produces  the  geometry  (BA) 
cards  read  by  CE-QUAL-ELV2 .  GIN  also  may  be  manipulated  to  change  the  com¬ 
puted  elevation-area-volume  table  to  more  closely  match  the  given  table. 

GIN  produces  an  active/inactive  cell  map  to  allow  screening  of  illegal  geo¬ 
metric  features. 

84.  Implicit  in  the  geometric  schematization  is  the  identification  of 
the  waterbody  centerline.  This  may  be  defined  as  the  dominant  flow  path, 
with  its  origin  at  the  upstream  boundary  and  it  may  follow  a  path  that 
coincides  with  the  thalweg.  All  longitudinal  distances  are  measured  along 
the  selected  centerline. 
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85.  The  choice  of  Ax  and  h  is  important  to  the  success  of  the 
simulation.  By  selecting  the  ratio  h/Ax  to  match  the  overall  estuary 
bottom  slope,  the  user  will  find  that  the  estuary  bottom  topography  may 
be  modeled  more  accurately.  While  relatively  small  values  of  Ax  and 
h  permit  finer  resolution,  the  resulting  large  number  of  grid  cells 
requires  more  computation  and  storage.  A  second  consideraion  in  the 
choice  of  Ax  and  h  is  the  basic  stability  criterion  that  At  ,  the 
integration  time  step,  be  smaller  than  the  volume  of  each  cell  divided 

by  the  flow  through  the  cell  (Equation  (8)).  This  criterion  also  affects 
economy  by  forcing  smaller  At's  as  cell  dimensions  are  decreased. 

86.  In  practice,  At  can  be  changed  for  a  particular  simulation, 
while  the  geometry  established  initially  remains  fixed.  Surprises  from 
an  implicit  choice  of  At  can  be  avoided  by  anticipating  the  limiting 
At  dictated  by  geometry.  In  most  cases,  this  limiting  At  will  occur 
at  tributary  or  upstream  inflows,  withdrawals,  or  downstream  outflows. 
Equation  (8)  can  be  evaluated  using  the  smallest  cell  volume  and  largest 
anticipated  flow  at  these  inflow  or  outflow  segments. 

87.  The  code  automatically  expands  the  grid  to  accommodate  rising 
and  falling  water  surfaces,  so  the  selected  grid  should  be  large  enough 
to  contain  the  largest  anticipated  volume.  The  grid  consists  of  active 
and  inactive  cells,  the  former  being  those  at  which  velocities,  tempera¬ 
tures,  salinities,  and  constituent  concentrations  are  computed;  and  the 
latter  being  those  at  which  no  computations  are  done.  Inactive  cells 
are  identified  by  inputing  their  widths  as  zero.  The  grid  must  satisfy 
the  following  rules: 

a.  It  must  be  at  least  two  active  cells  deep  at  every  segment. 

b.  It  must  be  at  least  three  active  cells  long  at  every  layer. 

c.  Steps  in  the  waterbody  bottom  profile  are  permitted  as  long 
as  each  continuous  layer  produced  is  at  least  three  cells 
long . 

d.  The  grid  must  include  a  layer  of  inactive  cells  at  the  top 
and  bottom  and  a  segment  of  inactive  cells  at  the  left  and 


right,  such  that  those  cells  representing  the  waterbody 
are  surrounded  on  all  sides  by  a  layer  or  segment  of  in¬ 
active  cells. 

e.  Cell  widths  must  not  decrease  from  bottom  to  top  in  any 
segment . 

Active  cells  are  those  which  may  contain  water,  or  may  not  at  various  times 
because  of  flood  events,  drawdown,  or  tidal  changes.  The  current  bound¬ 
aries  of  the  computation  on  the  grid  are  defined  by  the  parameters  IL 
and  IMAXM1  in  the  longitudinal  and  KT  and  KMAXM1  in  the  vertical. 

The  minimum  values  of  the  parameters  IL  and  KT  are  each  two,  indicating 
the  waterbody  is  at  its  maximum  volume.  As  drawdown  occurs,  KT  increases, 
so  that  the  water  surface  is  located  near  layer  KT  .  If  necessary,  CE- 
QUAL-ELV2  also  decreases  IL  so  that  segment  I  =  IL  is  always  at  least 
two  cells  deep. 

88.  Layers  are  identified  by  the  layer  number  K  ,  with  a  range  of  1  to 
KMAX  .  Each  of  these  has  an  elevation  associated  with  it,  so  that 

EL  =  h‘ (KMAX  -  KT  +  1)  -  Z  +  DTM  (43) 

is  the  water  surface  elevation,  where  DTM  is  the  distance  from  some 
reference  datum  to  the  bottom  of  the  grid.  CE-QUAL-ELV2  computes  water 
surface  elevations  in  terms  of  the  layer  number  K  and  the  local  elevation 
Z  ,  which  is  a  deviation  from  the  top  of  that  layer  rather  than  an  actual 
elevation.  Longitudinal  distances  are  given  in  terms  of  the  segment 
number  I  .  Mapping  these  segments  and  layers  back  onto  the  plan  and 
elevation  of  the  estuary  will  help  to  relate  CE-QUAL-ELV2  results  to  speci¬ 
fic  estuary  locations. 

89.  This  mapping  will  also  be  useful  in  establishing  the  locations 
of  inflows  and  outflows.  CE-QUAL-ELV2  uses  the  following  definitions: 


a.  Upstream  inflows  occur  only  at  I  *  IL  (left  boundary)  and 
are  defined  by  the  variables  QIN  ,  TIN  ,  and  CIN  . 

b.  Downstream  outflows  occur  only  at  I  =  IMAXM1  (right 
boundary)  and  are  defined  by  the  variable  QOUT  . 

c.  There  must  always  be  an  inflow  specified,  although  the  asso¬ 
ciated  flow  may  be  zero. 

d.  Tributaries  can  occur  in  any  segment  from  I  =  IL  to 
I  =  IMAXMl  ;  each  tributary  must  have  an  associated  flow, 
temperature,  salinity,  and  constituent  concentration. 

e.  Withdrawals  can  be  located  in  any  segment  and  layer,  and 
each  must  have  an  associated  flow. 

90.  The  estuary  orientation  may  be  specified  if  wind  speed  and  direc¬ 
tion  data  are  available  and  if  wind  stress  computations  are  desired.  The 
general  estuary  direction  is  determined  by  the  angle  the  centerline 
(x-axis,  positive  to  the  right  downstream)  makes  with  north. 

91.  As  presently  coded,  CE-QUAL-ELV2  will  handle  grid  sizes  up  to  40 
segments  long  and  50  layers  deep.  Larger  sizes  require  changes  in  dimension 
statements.  These  statements  are  marked  "DIMENS"  in  columns  73  to  80, 

and  instruction  i  for  changing  them  are  located  in  comment  statements  pro¬ 
ceeding  each  dimension  statement.  Dimension  statements  are  found  in 
the  main  program  and  in  subroutines  TRIDAG,  GRIDG,  GRID,  GRIDC,  and  GRIDX. 

Time-Varying  Boundary  Condition  Data 

92.  In  order  to  model  water,  heat,  and  salinity  budgets  accurately 
through  the  simulation  period,  it  is  necessary  to  consider  heat  ex¬ 
change  at  the  water  surface  and  mass,  heat,  and  salinity  advected  in  or 
out  of  the  waterbody.  Information  necessary  to  compute  these  fluxes  is 
provided  to  the  main  code  through  the  subroutine  TVDS,  which  is  an  acronym 
for  time-varying  data  selector.  This  subroutine  performs  two  tasks.  It 
is  called  once  in  order  to  read  the  time-varying  data.  On  each  succeed¬ 
ing  call,  TVDS  receives  from  the  main  code  the  current  simulation  time 
and  returns  to  the  main  code  the  value  of  eacn  of  the  parameters  specified 
in  its  call.  TVDS  must  be  modified  by  the  user  to  accommodate  his 
particular  data.  Minimum  data  requirements  are  flows,  temperatures, 


and  salinities  for  the  upstream  inflow  (QIN,  TIN,  SAIN)  and  each  tributary 
(QTRIB,  TTRIB,  STRIB) ;  flow  rates  for  each  withdrawal  (QWD) ;  and  the  sur¬ 
face  heat  exchange  parameters  CSHE  ,  SRO  ,  and  ET  .  If  evaporation  is 
to  be  computed,  the  dew  point  temperature  (TD)  is  required;  wind  speed 
(WA)  and  direction  (PHI)  are  needed  if  wind  stress  computations  are  to  be 
performed.  For  each  constituent  being  simulated,  inflow  (CIN)  and  tribu¬ 
tary  (CTRIB)  concentrations  are  also  required.  These  parameters  can  be 
supplied  at  any  time  interval  (including  irregular  intervals).  Daily  or 
more  frequent  data  is  recommended,  although  in  some  cases  steady-state 
simulations  made  with  constant  parameters  are  useful. 

93.  Right-hand  (open)  boundary  values  for  elevation  (ZR) ,  temperature 
(TR) ,  salinity  (SR),  and  constituent  concentration  (CR)  profiles  are  re¬ 
quired  and  are  also  provided  to  the  main  code  through  TVDS.  The  profile 
boundary  condition  data  may  be  observations  and  are  therefore  handled  in 
a  manner  similar  to  the  inflow  or  meteorological  records.  The  elevations 
may  also  be  observations  but  must  be  converted  to  the  system  used  in  the 
program;  i.e.,  the  elevations  must  be  in  terms  of  a  deviation  from  the 
current  top  layer  KT  .  Paragraphs  82  through  91  include  a  discussion  of 
the  elevation  convention,  and  Equation  43  gives  the  relationship  of  KT 
and  Z  to  an  elevation  relative  to  mean  sea  level. 

94.  The  right-hand  boundary  elevation  ZR  may  be  computed  more  con¬ 
veniently  from  the  tidal  amplitude  and  period  as: 

AZ  =  A  sin  (2Trt/t0)  +  ZI  (44) 

where  A  is  the  amplitude  (one  half  the  range),  t0  is  the  tidal  period 
(usually  12.45  hours),  and  ZI  is  the  initial  or  reference  elevation. 

95.  The  subroutine  TVDS  is  given  the  simulation  time  (variable  ELTM) 
which  is  measured  in  Julian  days  and  fractions  thereof.  ELTM  is  computed 
in  CE-QUAL-ELV2  as  TMSTRT  +  N  *  DLT/86400  ,  where  N  is  the  current  iter¬ 
ation  number.  Note  that  an  ELTM  of  1.75  means  January  1,  6  p.m.  and  that 
an  ELTM  of  72.0  is  midnight  March  12/March  13.  A  Julian  data  calendar  is 
given  in  Appendix  H. 


96.  In  the  most  general  case,  TVDS  has  available  from  its  first  call 
several  arrays,  each  one  of  which  contains  time-varying  dal  a  from  a  sep¬ 
arate  disk  or  tape  file.  In  the  example  for  the  Patuxent  Estuary,  one 
array  contains  all  the  flow  data  and  a  second  array  contains  the  right- 
hand  boundary  salinity  profile  data.  The  first  array  contains  daily  data, 
and  the  second  one  contains  data  taken  at  irregular  intervals.  These 
files  are  chronological,  which  is  essential  for  the  TVDS  algorithm.  Pairs 
of  times  associated  with  each  data  file  are  examined  sequentially  and  com¬ 
pared  to  ELTM  ,  beginning  with  the  previously  located  line  of  data.  When 
the  correct  line  is  located,  the  data  is  extracted  from  the  array  and  set 
equal  to  the  proper  subroutine  arguments.  In  the  example  shown,  both  data 
arrays  are  examined  in  this  manner  before  the  time-varying  boundary  condi¬ 
tion  data  are  returned  to  the  main  code.  The  subroutine  also  contains 
error-checking  routines. 

97.  It  is  important  that  the  prospective  user  study  the  TVDS  subroutine. 
Difference  in  data  availability  and  format  make  it  easier  for  the  user  to 
tailor  TVDS  to  his  data  than  to  attempt  to  create  a  generalized  TVDS  able 

to  handle  all  cases. 

Initial  Conditions,  Hydraulic  Parameters, 
and  Constituent  Reactions 

98.  Once  the  user  has  established  a  grid  to  represent  the  estuary  and 
assembled  time-varying  data  for  the  simulation  period,  initial  conditions 
to  start  each  simulation  must  be  selected.  These  consist  of  a  water  sur¬ 
face  elevation,  temperature,  and  salinity  corresponding  to  the  simulation 
starting  time.  This  information  is  coded  on  the  IC  input  card.*  If  the 
water  quality  constituent  computation  option  is  used,  initial  conditions 
for  each  constituent  are  required.  These  are  specified  on  the  CC  card. 

99.  The  initial  value  for  the  temperature  and  for  each  of  the  con¬ 
stituent  computations  is  a  single  number,  so  that  homogeneous  temperature 
and  constituent  fields  are  the  initial  conditions  at  which  the  simulation 
begins.  With  some  reprogramming,  temperature  and  constituent  fields  that 

*  Each  of  the  cards  or  card  images  used  as  input  to  CE-QUAL-ELV2  is  described 
in  Appendix  A  of  this  report. 
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vary  longitudinally  and  vertically  may  be  used.  This  option  is  already 
programmed  for  the  case  of  salinity.  A  single  value  may  be  used  to  ini¬ 
tialize  the  salinity  field  (IC  card),  or  a  separate  salinity  for  each 
active  cell  in  the  grid  may  be  input  (SI  cards) .  Appendix  A  gives  in¬ 
structions  for  selecting  these  options.  Although  the  simultaneous, 
implicit  solution  of  the  combined  momentum  and  continuity  equations 
establishes  a  consistent  flow  field  within  a  few  wave  travel  times  down 
the  estuary,  the  establishment  of  a  salinity  field  consistent  with  the 
applied  boundary  conditions  will  take  an  amount  of  time  equal  to  at 
least  a  residence  time.  For  this  reason,  it  is  recommended  that  a  lon¬ 
gitudinally  and  vertically  varying  salinity  field  be  used  for  the  initial 
condition.  There  are  two  ways  to  obtain  this  field:  (a)  using  observa¬ 
tions  or  (b)  computing  a  steady-state  field  and  using  the  resulting 
salinities  as  the  initial  condition  for  the  simulation  period  of  interest. 

100.  The  foregoing  comments  may  also  be  applied  to  the  initial  temper¬ 
ature  and  constituent  conditions.  However,  surface  heat  exchange  can 
reduce  substantially  the  amount  of  time  required  to  damp  out  the  initial 
condition  and  establish  a  temperature  field  consistent  with  the  boundary 
conditions.  The  same  may  be  said  of  any  constituents  that  have  a  decay 
term  on  the  right-hand  side  of  the  transport  equation. 

101.  Hydraulic  parameters  consist  of  ,  Az  »  and  the  Chezy 
coefficient.  These  parameters  are  defined  in  Part  II  and  discussed  in 
Part  III  of  this  guide,  and  suggested  values  are  coded  in  DATA  statements. 

and  A^  may  be  varied  by  several  orders  of  magnitude,  with  the 

following  exception.  A  ,  the  vertical  dispersion  of  momentum  coefficient, 

—  6  2 

has  the  limits  of  the  molecular  value  (1.5  x  10  itT/s)  as  a  minimum  and  a 

value  related  to  the  horizontal  layer  thickness  and  the  time  step  as  a 

maximum  (A^  <  h2/2At).  The  latter  limit  is  derived  from  the  computation 

time  step  limit,  Equation  (8).  The  value  of  t\?  varies  spatially  and 

temporally  within  these  limits  and  is  computed  by  the  code  from  the 

Richardson  number,  which  is  the  local  ratio  of  buoyancy  to  shear.  The 

value  of  A^  used  in  the  computation  is  a  function  of  a  base  value  for 

this  parameter  A  and  the  Richardson  number,  so  that 
zo 


49 


(45) 


A  =  A  (1  +  lORi)  ** 
z  zo 

where  Ri  is  the  Richardson  number.  It  is  necessary  to  set  only  A  , 

zo 

as  the  code  automatically  keeps  the  value  of  A^  within  the  molecular 
and  computational  limits.  A^  is  presently  set  at  a  single  value  that 
is  constant  over  time  and  space,  but  the  numerical  formulation  allows  it 
to  vary  spatially  and  temporally  with  additional  coding. 

102.  The  Chezy  coefficient  is  also  space-  and  time-invariant,  but  it 
can  be  made  variable  with  some  reprogramming.  Boundary  friction  varies 
spatially  with  the  amount  of  side  and  bottom  area  exposed,  as  well  as 
temporally  with  the  flow  field.  Chezy  may  be  varied  by  a  factor  of  four 
from  the  value  given  in  the  CE-QUAL-ELV2  DATA  statement. 

103.  The  temperature  equation  dispersion  parameters  may  be  varied. 

These  parameters  are  D  and  D  ,  the  longitudinal  and  vertical  disper- 

X  z 

sion  coefficients.  Like  A^  ,  is  temporally  and  spatially  invariant, 

and  its  value  is  controlled  by  the  FORTRAN  variable  DXI  .  D  varies 

z 

over  time  and  space  and  is  also  computed  from  the  Richardson  number  and  a 
base  value  Dzq  .  Like  the  momentum  equation  parameter  Az  ,  Dz  oper¬ 
ates  within  a  range  from  a  minimum  molecular  value  to  a  computational 
maximum  of  h2/2At  . 

104.  The  distribution  of  solar  radiation  in  the  vertical  is  controlled 
by  the  parameters  8  and  Y  (FORTRAN  BETA  and  GAMMA)  as  follows: 

Hg(z)  =  (1  -  8)  Hs(z  =  0)“YZ  (46) 

where 

Hg  solar  radiation  rate  at  a  depth  z  (°C  m3,s  Vm  ) 

8  fraction  of  H  absorbed  at  the  water  surface 

s 

Y  attenuation  rate 

Note  that  these  parameters  control  only  the  distribution  of  solar  radiation 
in  the  vertical  and  that  overall  surface  heat  exchange,  including  solar 
radiation,  is  summarized  in  the  parameters  CSHE  and  ET  .  If  8  =  3  » 
all  the  solar  radiation  is  absorbed  in  the  current  top  layer  (K  =  KT)  , 
and  the  value  of  Y  is  arbitrary. 


Constituent  Reactions 


105.  For  each  constituent  specified,  constituent  reaction  rates  and 
internal  sources  and  sinks  need  to  be  coded.  For  the  user  who  desires 
only  hydrodynamic  and  temperature  simulations,  this  section  of  the  user 
guide  may  be  skipped.  For  a  detailed  description  of  the  constituent 
transport  computations,  the  user  may  refer  to  LARM2  Developments  1979-1980 
(Edinger  and  Buchak  1980b). 

106.  The  variable  RR(JC,M)  is  the  rate  at  which  constituent  JC  is 
transferred  to  constituent  M  and  may  be  dependent  on  a  number  of  spatially 
and  temporally  varying  variables,  such  as  temperature.  For  the  case  of 

M  =  JC  ,  RR  is  negative  and  is  the  linear  decay  rate.  RR  is  expressed 
in  units  of  s-1  . 

107.  The  variable  HN  (I,K)  is  the  spatially  and  temporally  varying 

source  and  sink  term  for  constituent  JC  .  HN  (I,K)  has  the  units 

—  1  3  —  1 

mg  1  -m  -s  .  The  user  is  required  to  augment  this  array  by  the  amount 
of  material  added  to  or  subtracted  from  constituent  JC  due  to  the  reac¬ 
tions  specified  by  RR  . 

108.  The  user  must  examine  DO  loop  1840  in  the  CE-QUAL-ELV2  code  and 
modify  it  to  compute  RR  and  HN  as  required  for  his  constituents.  The 
user  may  also  specify  different  sources  and  sinks  in  the  surface  or  bottom 
layers  by  separate  statements  for  layers  K  =  KT  or  K  =  KB  .  The  compu¬ 
tation  for  each  constituent  concentration  (DO  loop  1820)  begins  with  an 
initialization  of  the  HN  array.  Then,  for  each  location  I,K  the  reac¬ 
tion  rate  (RR)  array  is  computed,  one  FORTRAN  statement  for  each  constit¬ 
uent  JC  is  transferred  to  or  from  constituent  M  .  The  remainder  of  the 
computation  requires  no  user  changes  and  includes  the  augmenting  of  HN 

by  external  sources  and  sinks,  the  retrieving  of  the  transport  coefficient 
vectors,  and  the  call  to  the  subroutine  TRIDAG  for  layer-by-layer  solution 
of  the  transport  equation. 

Input,  Output,  and  Computer 
Resource  Requirements 

109.  The  input  data  stream  is  described  in  Appendix  A,  and  an  example 
data  set  for  the  Patuxent  Estuary  is  given  in  Appendix  B.  CE-QUAL-ELV2 
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output  for  this  data  set  is  given  in  Appendix  C.  The  output  shown  there 
is  for  IFORM  =  0  ,  that  is,  compressed  and  shortened  for  80-character- 
width  printers.  The  normal  form  of  the  output  is  for  132-character-width 
printers.  Although  an  example  of  this  more  complete  output  is  not  given, 
the  user  may  quickly  make  his  own  by  changing  the  IFORM  parameter  (PR 
card,  field  1)  from  0  to  1.  T  e  additional  information  generated  by  the 
longer  output  is  as  follows: 

a.  Three  geometry  tables:  a  grid  showing  active  and  inactive 
cells,  a  grid  showing  cell  widths  as  augmented  by  CE-QUAL-ELV2 
(widths  assigned  to  inactive  cells) ,  and  a  elevation-area- 
volume  table. 

b.  The  initial  conditions  for  Z2  ,  U  ,  W  ,  and  S2  (and 
C2  ):  the  velocity  fields  are  always  initialized  to  zero 
at  the  start  of  the  computation,  Z2  and  T2  are  ini¬ 
tialized  through  the  IC  card,  C2  through  the  CC  card, 
and  S2  through  the  IC  or  SI  cards. 

c.  Eight  additional  segment  results  for  Z2  ,  V  ,  W  ,  and 
S2  (and  C2  ) ,  using  the  additional  columns  available  on 
11-  by  14-inch  paper. 

110.  A  second  type  of  output  is  available  from  the  CE-QUAL-ELV2  code,  and 
it  is  useful  in  postprocessing  applications  such  as  the  vector  plotting  code 
W  .  Cell  coordinate  positions  are  written  unformatted  to  TAPE61,  followed 
by  Z2  ,  U  ,  and  W  fields  at  times  selected  on  the  PL  card  (see 
Appendix  A).  This  information  is  used  by  W  to  plot  displacement  vectors 
on  a  scaled  grid.  However,  by  saving  additional  information  on  TAPE61, 

time  series  plots  of  other  variables  may  be  obtained,  as  well  as  statisti¬ 
cal  analyses  of,  for  example,  outflow  temperatures  or  velocity  records  at 
a  particular  estuary  location. 

111.  The  code  and  test  data  as  supplied  requires  320,000  octal  words  of 
storage  on  a  Controled  Data  Corporation  175.  This  storage  requirement  de¬ 
pends  on  both  the  grid  size  and  the  amount  of  time-varying  boundary  condi¬ 
tion  data  stored.  CE-QUAL-ELV2  requires  0.00023s  of  cpu*  time  for  each  time 


*  central  processor  unit. 
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step  iteration  for  each  active  cell  at  a  high  level  of  compiler  optimi¬ 
zation.  The  Patuxent  Estuary  example  uses  124.72s  of  cpu  time  for  its 
452  active  cells  for  simulation  of  1200  steps,  plus  11.5s  for  compilation. 

Patuxent  Estuary  Example 

112.  The  problem  here  was  to  determine  how  to  distribute  heat  from 
the  Chalk  Point  Power  Plant  within  the  longitudinal-vertical  structure  of 
the  Patuxent  Estuary.  The  characteristics  of  the  estuary  and  the  ratio¬ 
nale  for  applying  CE-QUAL-ELV2  to  this  particular  problem  are  presented  in 
Part  III  of  this  report.  This  application  presented  the  opportunity  to 
verify  the  computations  with  observed  velocities,  temperature,  and 
salinities . 

113.  This  section  of  the  report  describes  the  steps  in  applying 
CE-QUAL-ELV2 ,  with  particular  emphasis  on  the  coding  involved.  This  appli¬ 
cation  is  typical  because  some  features  of  this  problem  could  not  be  handled 
by  the  general  boundary  condition  routines  of  CE-QUAL-ELV2 .  TVDS  was  re¬ 
written  to  supply  the  particular  time-varying  boundary  condition  data  to  the 
main  code;  however,  the-  user  will  find  elements  of  many  applications  in  this 
Patuxent  Estuary  example  and  is  encouraged  to  study  it. 

114.  Figure  3  shows  the  Patuxent  Estuary.  Cronin  and  Pritchard  (1975) 
present  dimensional  data  for  the  Patuxent,  including  average  widths  for 
each  nautical  mile  in  the  longitudinal  direction  and  for  each  metre  in 
the  vertical.  These  data  were  used  directly  in  creating  a  finite  differ¬ 
ence  representation  of  the  Patuxent,  with  a  Ax  of  one  nautical  mile 
(1852  m)  and  an  h  of  one  metre;  therefore  the  preprocessors  GEDA  and 
GIN  were  not  needed  in  this  application.  The  resulting  grid  with 
IMAX  =  36  and  KMAX  =  31  is  shown  in  Figure  4.  The  segments  are  also 
mapped  onto  the  Patuxent  Estuary  plan  in  Figure  5. 

115.  The  level  of  zero  depth  in  the  tables  given  by  Cronin  and 
Pritchard  is  mean  low  water.  The  parameter  DTM  was  assigned  a  value 
of  -30  m  so  that  the  computation  of  EL  (Equation  (43))  yields  0  m  when 
KT  is  2  and  Z  is  0  m.  EL  ,  therefore,  becomes  the  deviation  from 
mean  low  water  at  the  open-water  (Chesapeake  Bay)  boundary. 
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116.  Following  the  establishment  of  the  finite  difference  grid,  in¬ 
flows,  tributaries,  withdrawals,  and  the  open-water  boundary  were  identi¬ 
fied  and  located  on  the  grid.  These  are  as  follows: 

a.  Patuxent  River:  QIN  ,  TIN  ,  SAIN  — upstream  inflow 

b.  Chalk  Point  Power  Plant  intake:  QWD(l)  — segment  13, 
layer  2 

c.  Chalk  Point  Power  Plant  discharge:  QTRIB(l)  ,  TTRIB(l)  , 

STRIB(l)  —segment  11 

c.  Chesapeake  Bay:  ZR  ,  TR  ,  SR  — open-water  boundary. 

The  discharge  temperature  and  salinity  are  dependent  on  the  temperature 
and  salinity  adjacent  to  the  intake,  as  well  as  on  the  amount  of  heat 
added  by  condenser  cooling.  This  dependency  can  be  handled  in  a  number 
of  ways;  for  example,  the  temperature  and  salinity  at  the  Chalk  Point 
Power  Plant  intake  could  be  passed  to  TVDS  through  the  call  parameters, 
with  TTRIB(l)  and  STRIB(l)  set  in  TVDS.  A  simpler  method  was  adopted. 
TVDS  supplied  the  main  code  with  the  plant  temperature  rise,  and  the  fol¬ 
lowing  two  statements  were  added  to  the  main  code  following  the  call  to 
TVDS: 


r 

TTRIB(l)  =  TTRIB(l)  +  T2  IWD(l),  KWD(l) 

(47) 

STRIB(l)  =  S2  IWD(l),  KWD(l) 

(48) 

The  variable  TTRIB(l)  on  the  right-hand  side  of  the  first  statement  is  the 
plant  temperature  rise  assigned  in  TVDS.  This  code  is  the  only  modification 
required  in  CE-QUAL-ELV2  for  the  Patuxent  Estuary  application,  other  than  the 
rewriting  of  TVDS. 

117.  The  time-varying  data  consist  of  two  tapes:  observed  inflows  and 
boundary  salinities  and  generated  boundary  tide.  The  observations  of  daily 
QIN  values  were  recorded  for  July  and  August  1978,  and  five  right-hand 
boundary  salinity  profiles  were  made  for  a  period  beginning  17  July  and 
ending  18  September.  The  boundary  tide  was  synthesized  as  follows  using 
a  period  of  12. A5  hours  and  an  extrapolated  amplitude  for  Solomons,  Md . : 


ZR  =  -0.24  *  SIN(12. 1122  *  ELTM)  +  ZI 


(49) 


where  ZI  is  the  initial  water  surface  elevation  relative  to  layer  two. 

118.  All  other  boundary  condition  data  were  assumed  to  be  constant. 
Note  that  this  simulation  was  designed  to  obtain  excess  temperatures  above 
a  base  of  20°C.  These  were  obtained  by  setting  the  inflow  and  open-water 
boundary  temperatures  to  20°C  as  well  as  the  equilibrium  temperature  ET 
to  20°C.  The  coefficient  of  surface  heat  exchange  CSHE  was  set  to  a 
nominal  value  of  120  Btu/(ft2  *day°F) . 

119.  Initial  conditions  for  the  simulation  correspond  to  the  values 
of  the  boundary  condition  data  for  the  selected  starting  date,  Julian 
day  205  (24  July  1978).  Initial  salinities  were  taken  from  the  same  set 
of  observations  that  produced  the  boundary  salinities.  These  data  per¬ 
mitted  the  salinity  field  to  be  initialized  with  vertical  and  longitudinal 
gradients.  The  initial  temperature  was  set  to  20°C  to  be  consistent  with 
the  excess  temperature  computations  on  the  20°C  base. 

120.  An  initial  series  of  simulations  was  made  to  check  all  the 
boundary  condition  data  as  supplied  to  the  main  code  by  TVDS  and  to  help 
select  a  At  .  An  acceptable  value  for  At  proved  to  be  747s,  which  was 
large  enough  for  economical  simulations,  yet  small  enough  to  resolve  the 
tidal  cycle  that  drives  the  computation  (747s  is  1/60  of  the  tidal  period 
of  12.45  hours).  The  final  simulation  included  a  long  period  of  initial¬ 
ization  (1080  time  steps,  or  9.3  days)  followed  by  two  entire  tidal 
cycles  (120  additional  steps),  during  which  results  were  printed  every 
five  time  steps. 
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APPENDIX  A:  INPUT  DATA  DESCRIPTION 
FOR  CE-QUAI.-ELV2 


General 

Description  for  each  of  the  cards  (or  card  images)  used  as  input  to  CE-QUAL- 
ELV2  are  presented  on  the  following  pages.  There  are  three  data  types: 


alphanumeric 

real  (may  be  left-  or  right-justified) 
integer  (must  be  right-justified) 

There  are  eleven  data  fields  per  card: 

Columns 


1-2 

3-8 

9-16 

17-24 

25-32 

33-40 

41-48 

49-56 

57-64 

65-72 

73-80 


Field 

0 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Length 

2 

6 

8 

8 

8 

8 

8 

8 

8 

8 

8 


This  format  is  similar  to  that  used  in  the  Hydrologic  Engineering  Center 
codes . 


A! 


Title  Card  (Required) 


Title  cards  may  be  used  to  identify  applications,  simulations,  or  other 
parameters  as  required.  These  cards  do  not  otherwise  affect  the 
computation. 


CARD  T 1 


The  text  input  with  this  card  appears  in  the  header  with  each  velocity, 
temperature,  salinity,  and  constituent  field  snapshot,  as  well  as  in  the 
overall  header  printed  once  for  each  simulation. 


Field 

Parameter 

Value 

Description 

0 

AID 

T1 

alphanumeric: 

card  identification 

1-9 

TITLE( 1 , J)  , 
J=L , 18 

alphanumeric : 
tion  text 

simulation  identifies 

CARD  T2 


The  text  input  using  this  card  always  appears  in  the  overall  simulation 
header  and  appears  in  the  first  snapshot  header  if  IFORM  *  1  . 


Field 

Parameter 

Value 

Description 

0 

AID 

T2 

alphanumeric: 

card  identification 

1-9 

TITLE(2 , J) , 
J=1 , 18 

alphanumeric: 
tion  text 

simulation  identifica- 

CARD  T3 


The  text  input  using  this  card  always  appears  in  the  overall  simulation 
header  and  appears  in  the  first  snapshot  header  if  IFORM  =  1  . 


Field 

Parameter 

Value 

Description 

0 

AID 

T3 

alphanumeric : 

card  identification 

1-9 

TITLE(3, J) , 

alphanumeric : 

simulation  identifica 

J=1 , 18 

tion  text 

CARD  T4 


The  text  input  using 

this  card 

always  appears  in 

the  overall  simulation 

header 

and  appears  in 

the  first 

snapshot  header 

if  IFORM  =  1  . 

Field 

Parameter 

Value 

Descr i ption 

0 

AID 

T4 

alphanumeric : 

card  identification 

1-9 

TITLE(4 , J) 

* 

alphanumeric : 

simulation  identifica- 

J=l, 18 

tion  text 

Note:  The  variable  TITLE  is  dimensioned  to  accommodate  18  four-charac¬ 
ter  words  per  line.  This  configuration  will  allow  CE-QUAL-ELV2  to 
run  on  machines  that  have  at  least  four  characters  per  word,  such  as 
IBM  computers. 


Geometry  Card  (Required) 


CARD  GE 


The  geometry  card  specifies  the  grid  size  (number  of  segments  and  layers), 
the  length  and  height  of  the  individual  cells,  and  the  location  of  the 
grid  relative  to  a  reference  datum. 


Field 

Parameter 

Value 

Description 

0 

AID 

GE 

alphanumeric:  card  identification 

1 

IMAX 

>5 

integer:  number  of  segments 

2 

KMAX 

>4 

integer:  number  of  layers 

3 

DLX 

+ 

real:  Ax  ,  cell  length,  m 

4 

HIN (1,1) 

+ 

real:  h  ,  cell  height,  m 

5 

DTM 

-,0,+ 

real:  distance  from  reference  datum 

to  lower  grid  boundary,  m 

Note:  CE-QUAL-ELV2  requires  a  grid  with  IMAX>5  and  KMAX>4  .  A  typi¬ 
cal  grid  has  IMAX  in  the  10-to-36  range,  KMAX  in  the  15-to-31 
range.  Grid  sizes  larger  than  36  *  31  require  changing  the 
code's  dimension  statements.  Instructions  for  doing  so  are 
given  in  comment  cards  preceding  each  dimension  statement, 
each  of  which  is  marked  "DIMENS"  in  columns  73-80. 

Typical  cell  sizes  are  DLX  in  the  0.5-to-10-km  range  and 
HIN(1,1)  in  the  0.5-to-30-m  range.  HIN(I,i)  is  read  as  the 
cell  height  for  cell  1=1,  K  =  1  ,  which  is  then  applied 
to  all  other  cells.  DLX  likewise  is  the  length  of  each 
cell . 


A4 


Time  Card  (Required) 

CARD  TM 

This  card  determines  the  simulation  length  by  specifying  the  temporal 


iteration 

size  and  the 

number  of 

iterations. 

Field 

Parameter 

Value 

Description 

0 

AID 

TM 

alphanumeric:  card  identification 

1 

DLT 

+ 

real:  At  computation  time  step,  s 

2 

NSTEPS 

+ 

integer:  number  of  iterations 

Note:  The  simulation  length  is  NSTEPS  *  DLT  ,  in  seconds.  The  value 

of  DLT  may  be  determined  approximately  from  Equation  8.  However 
some  experimentation  may  be  necessary  to  determine  a  value  of  DLT 
at  which  the  computation  is  stable.  Also,  to  accurately  model  a 
tidal  boundary  condition  with  a  period  of  12.45  hours,  a  DLT  of 
no  more  than  15  minutes  should  be  used. 


CARD  PL 


This  card  specifies  whether,  and  at  what  frequency,  elevation  and  velocity 
fields  are  saved  for  postprocessing  by  WE,  the  vector  plotting  code. 


Field 

Parameter 

Value 

Description 

0 

AID 

PL 

alphanumeric:  card  identification 

1 

I  PLOT 

0 

integer:  no  postprocessing  information 
written  to  TAPE61 

1 

integer:  postprocessing  information 
written  to  TAPE61 

2 

Ml 

0,+ 

integer:  information  is  written  be¬ 

ginning  at  iteration  Ml  ;  if  IPLOT  =  0 
this  parameter  has  no  effect 

3 

M2 

0,+ 

integer:  information  is  written  every 

M2  iterations;  if  IPLOT  =  0  ,  this 
parameter  has  no  effect 

Note : 

If  IPLOT  =  1  ,  Ml  =  960  ,  and  M2  =  96  ,  for  example,  the  eleva¬ 
tion  and  velocity  fields  are  written  to  TAPE61  every  96  iterations, 
beginning  at  iteration  960. 

Print  Card  (Required) 


CARD  PR 


This  card  specifies  the  form  of  the  printed  output,  as  well  as  the  - 
quency  of  the  velocity,  temperature,  salinity,  and  constituent  field 
snapshots . 


Field 

Parameter 

Value 

0 

AID 

PR 

1 

IFORM 

0 

1 

2  Nl  0,+ 

3  N2  0,+ 

4  ID  LAC,  0 

1 


Description 

alphanumeric:  card  identif icatir n 

integer:  output  shortened  (cell 

widths  and  volume-area-elevation 
table  suppressed)  and  compressed 
to  fit  an  8%-inch  page  width;  snap¬ 
shot  output  is  for  9  segments 
specified  on  10  card 

integer:  normal  output  using  en¬ 

tire  132-column  page  width  (14 
inches) ;  snapshot  output  is  for 
17  segments  specified  on  II  and 
12  cards 

integer:  snapshots  are  printed 

beginning  at  iteration  Nl 

integer:  snapshots  are  printed 

every  N2  iterations 

integer:  the  printing  of  the  in.;  • 

vidua  1.  terms  in  the  momentum  equa¬ 
tion  is  suppressed 

integer:  the  printing  of  the  ir.': 

vidual  terms  in  the  momentum  equa¬ 
tion  occurs  according  to  the  co: 
parameters  IFORM,  Nl  ,  N2  , 

II  ,  and  12 


-• 


Note:  The  parameters  Nl  and  N2  are  used  like  the  parameters 

and  M2  on  the  PL  card.  The  segments  for  which  elevati  'n 
velocities,  temperatures,  salinities,  and  constituents  art 
printed  in  the  snapshots  are  specified  on  the  10  .  11 

12  cards.  Very  few  simulations  will  run  with  ID  I  AC  -  1 


A  7 


These  cards  permit  the  user  to  select  which  segment  results  are  to  be 
printed  for  IFORM  =  0  and  IFORM  =  1  .  Nine  segments  may  be  selected 
for  IFORM  *  0  (short  form)  and  17  segments  for  IFORM  =  1  (long 
form) . 

CARD  10 

Field  Parameter 

0  AID 

1-9  10(1), 

1=1,9 

CARD  II 

Field  Parameter  Value 

0  AID  II 

1-10  11(1),  1<I1<IMAX 

1=1,10 

CARD  12 

Field  Parameter  Value  Description 

0  MD  12  alphanumeric:  card  identification 

1-7  11(1),  1<I1<IMAX  integer:  segment  number 

1=11,17 


Value  Description 

10  alphanumeric:  card  identification 

1<I0<IMAX  integer:  segment  number 


Note:  If  IMAX <  9  or  IMAX <  17  ,  segment  numbers  may  be  repeated  to 

complete  the  card  requirements. 


• 

- 

P.  ' 

Initial  Condition  Card  (Required) 

CARD 

IC 

:c 

P- ... 

This 

card  defines  the  starting  time  in  Julian  days  and  the  corresponding 

• '  .  •  ''  m  , 

initial  water  surface  elevation  and  temperature. 

•  .  t  ' 

Field 

Parameter  Value 

Description 

I 

• 

*•  4 

0 

AID 

IC 

alphanumeric:  card  identification 

1 

TMSTRI 

0<TMSTRT<366 

real:  simulation  starting  time. 

Julian  days 

| 

2 

KT 

KMAX-2<KT<2 

integer:  layer  number  for  initial 

4 

• 

water  surface  elevation 

3 

ZI 

-0.8*HIN(I, 1)<ZI 

real:  initial  water  surface  elevation 

<0.25*HIN(1 , 1) 

relative  to  layer  KT  ,  m  (positive 

down) 

« 

r 

4 

TI 

0<TK100 

real:  initial  temperature,  C 

5 

SI 

-1 

real:  initial  salinity  is  specified  on 

, 

SI  cards  which  follow  BA  cards 

iWM 

0<SI<40 

real:  initial  salinity  of  SI  applied 

-• 

to  all  active  cells,  ppt 

fl 

Note: 

KT  and 

ZI  give  the  water  surface  elevation,  TI  the  temperature, 

and  SI 

the  salinity  for 

corresponding  to  the  time  TMSTRT  .  Values 

of  3  for 

KT  and  0.1  m  for  ZI  ,  for  example,  place  the  initial 

water  surface  elevation  0. 

1  m  below  the  top  layer  3.  The  initial 

temperature  field  is  isothermal  at  TI  and  isohaline  at  SI  unless 

SI  =  -1 

The  simulation 

begins  at  TMSTRT  and  ends  at  TMSTRT  + 

4 

NSTEPS  * 

DLT/86400  ,  days. 

• 

< 

• 

-  •¥ 

1 

• 

• 

1 

CARD  CD 


uent  Definition  Card  (Required) 


This  card  specifies  the  number  of  water  quality  constituents  to  be  simu¬ 
lated  and  also  permits  these  transport  computations  to  be  turned  off  for 
a  particular  simulation. 


Field  Parameter  Value  Desc.r  i  ption 


0 

1 


2 


AID 


NC 


ICC 


CD  alphanumeric:  card  identification 

0<NC<10  integer:  number  of  constituents  to  be 

simulated;  if  NC >  0  ,  the  following 
card.  Constituent  Initial  Concentration 
Card,  is  required 

0  integer:  constituent  transport  compu¬ 

tations  suppressed 


1 


integer:  constituent  transport  compu¬ 

tations  performed 


Note:  The  user  generally  will  complete  a  hydrodynamic,  temperature,  and 
salinity  simulation  to  his  satisfaction  (ICC  =  0)  before  going  on 
to  the  water  quality  constituent  computations  (ICC  =  1). 

Water  quality  reaction-iteration  rates  also  must  be  specified  for 
each  water  quality  constituent.  Since  these  rates  are  usually  a 
function  of  temperature,  they  must  be  computed  in  CE-QUAL-ELV2 .  The 
location  for  the  user  to  insert  this  computation  is  discussed  in 
paragraph  108  of  this  guide. 


J> 


AID 


CARD  CC 


This  card  contains  the  initial  concentration  for  each  of  the  NC  con¬ 
stituents  and  is  required  if  NC  >  0  . 


Field 

Parameter 

Value 

Description 

0 

AID 

CC 

alphanumeric : 

card  identification 

1-10 

CI(JC) 

JC=1 ,NC 

CI(JC)>0 

real:  initial 
each  of  the  NC 

concentrations  for 
constituents 

Note:  CE-QUAL-ELV2  initializes  the  entire  estuary  to  Cl(JC) 


Meteorological  Parameter  Card  (Required) 


CARD  MP 


This  card  specifies  the  estuary  orientation  for  wind  stress  computations 
and  also  whether  evaporation  is  included  in  the  water  budget  computations 


Field 

Parameter 

Value 

Description 

0 

AID 

MP 

alphanumeric:  card  identification 

1 

PHIO 

0<PHIO<2tt 

real:  the  angle  between  the  positive 

x-axis  (estuary  centerline  in  the 
flow  direction)  and  north,  radians 

2 

IEVAP 

0 

integer:  evaporation  rates  not  com¬ 
puted  and  not  included  in  water  bud¬ 
get  computations 

1 

integer:  evaporation  rates  computed 

and  included  in  water  budget  computa¬ 
tions 

Note:  For  an  estuary  flowing  from  east  to  west,  PHIO  =  it / 2  ;  the  south¬ 

west  to  northeast,  PHIO  =  5tt/4  .  The  convention  is  similar  to 
the  meteorologist's  convention  for  wind  direction. 

Evaporation  rates  are  sometimes  given  in  the  inflow  record  as 
QINnet  =  QIN-QET  .  This  is  the  case  for  which  IEVAP  should  be 
set  to  0.  The  contribution  of  evaporation  to  surface  heat  exchang 
is  always  considered,  regardless  of  the  value  of  IEVAP  . 


CARD  TR 


This  card  specifies 
Field  Parameter 

0  AID 


the  number  of  estuary  tributaries. 

Value  Description 

TR  alphanumeric:  card  identification 


1  NTRIB  0<NTRIB<10  integer:  number  of  estuary  tribu¬ 

taries;  if  NTRIB >  0  ,  then  the 
following  card.  Tributary  Location 
Card,  is  required 


Note:  For  each  tributary,  an  inflow  and  inflow  temperature  and  salinity 

need  to  be  supplied  by  the  subroutine  TVDS. 


Tributary  Location  Card  (Optional) 


CARD  TL 


This  card  gives  the  estuary  segment  numbers  at  which  tributary  inflows 
enter  and  is  required  if  NTRIB  >  0  . 


Field 


Parameter 


Value 


0 


AID 


TL 


Description 

alphanumeric:  card  identification 


1-10  ITRIB(J),  2<ITRIB(J)  integer:  segment  number  at  which  each 
J=l, NTRIB  <IMAX-1  tributary  enters 


Note:  Tributary  flows  are  entered  in  the  segment  specified  here  and  at 
a  layer  that  corresponds  to  their  density,  which  is  computed  by 
CE-QUAL-ELV2  and  which  is  shown  on  the  output  as  KTRIB  .  Trib¬ 
utaries  above  the  current  upstream  boundary  segment  IL  on  any 
particular  time  step  are  combined  with  the  upstream  inflow. 


Withdrawal  Card  (Required) 


CARD  WD 


This  card  specifies  the  number  of  withdrawals  from  the  estuary. 


Field 

Parameter 

Value 

Description 

'  - - “  -  J  ' 

0 

AID 

WD 

alphanumeric:  card 

identification 

1 

NWD 

0<NWD< 10 

integer:  number  of 
drawals;  if  NWD  >  0 
lowing  cards,  WI  and 
required 

estuary  with- 
,  then  the  fol 
WK,  are 

Note:  For  each  withdrawal,  an  outflow  rate  needs  to  be  supplied  by 

subroutine  TVDS. 


CARD  WI 


The  card  specifies  the  longitudinal  location  (segment)  for  each  with 
drawal  and  is  required  if  NWD  >  0  . 

Field  Parameter  Value  Description 

0  AID  WI  alphanumeric:  card  identification 

1—10  IWD(J),  2<IWD(J)  integer:  segment  number  for  each 

J=1,NWD  <IMAX-1  of  NWD  withdrawals 


Note:  The  layer  number  for  each  withdrawal  is  specified  on  the  WK 

card,  which  follows. 


CARD  WK 


This  card  specifies  the  vertical  location  (layer)  for  each  withdrawal 
and  is  required  if  NWD  >  0  . 

Field  Parameter  Value  Description 

0  AID  WD  alphanumeric:  card  identification 

1-10  KWD(J) ,  2<KWD(J)  integer:  layer  number  for  each  of 

J= 1 ,NWD  <KMAX-1  NWD  withdrawals 


Note:  The  segment  number  for  each  withdrawal  is  specified  on  the  pre¬ 
ceding  card.  Each  individual  withdrawal  must  have  its  segment 
and  layer  specified  in  the  same  field  on  the  WI  and  WK  cards, 
respectively. 

A  single  withdrawal  spanning  a  number  of  layers  may  be  speci¬ 
fied  by  breaking  it  into  several.  Withdrawals  above  the  cur¬ 
rent  upstream  boundary  segment  IL  or  above  the  current  water 
surface  layer  KT  are  ignored. 


Bathymetry  Cards  (Required) 


CARDS  BA 


These  cards  contain  the  cell  widths  for  each  of  the  IMAX *  KMAX  cells, 
which  are  taken  from  cross-section  drawings  and  represent  an  average 
value  for  each  cell.  GEDA  produces  widths  as  a  function  of  elevation. 
GIN  (see  Appendix  F)  uses  GEDA  output  and  produces  BA  cards  used  by 
CE-QUAL-ELV2 . 


F  ield 

Parameter 

Value 

Description 

0 

AID 

BA 

alphanumeric: 

card  identification 

1-10 

B  (I ,  K) 

B(I,K)>0 

real:  cell  widths,  m  (left  to  right, 

top  to  bottom,  filling  each  BA  card) 

Note:  For  the  grid  boundary  segments  (1=1  and  I  =  IMAX)  and  boundary 

layers  (K  =  1  and  K  =  KMAX),  B  must  be  set  to  zero;  inactive 
cells  also  are  identified  by  setting  their  widths  to  zero. 

The  read  statement  for  these  cards  is 
READ(5 , 550)  ((B(I,K) , K=1 , KMAX) , 1=1 , IMAX) . 


Salinity  Cards  (Optional) 


CARDS  SI 


These  cards  specify  the  initial  salinity  field  and  are  required  if  the 
SI  parameter  on  the  IC  card  is  -1. 


Field 


Parameter 


AID 


Value 


SI 


Description 

aiphnanumeric :  card  identification 


1-10  S1(I,K)  0<S1(I,K)<40  real:  initial  salinities  for  each 

active  cell  (0  for  inactive  cells), 
ppt  (left  to  right,  top  to  bottom, 
filling  each  SI  card) 


Note:  The  read  statement  for  these  cards  is 

IF (SI . EQ . - 1 . )  READ ( 5 , 550)  ((SI (I ,K) ,K=  1 ,KMAX) ,1=1, IMAX) . 


A 19 


TVDS  Data  Cards  (Optional) 


These  cards  contain  the  time-varying  boundary  condition  data  for  the 
simulation.  Since  they  are  read  by  the  user-written  subroutine  TVDS, 
no  format  is  specified  here.  A  discussion  of  the  time-varying  boundary 
condition  data  and  TVDS  can  be  found  in  paragraphs  92-97. 
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APPENDIX  C:  EXAMPLE  CE-QUAL-ELV2  OUTPUT  DATA 

L  A  E  N  --  LATERALLY  AVEKAGEO  ESTUARY  MODEL 

VERSION  TkG  --  OCTOBER  Hbl  —  aAYnE.  PENNSYLVANIA 

USER  GolUt  Example  --  PATLXENT  ESTUAPy.  MARYLAND  --  13/1A/M 
UPSTREAM  INFlG..  FOwER  PLANT  INTAKE  AN 0  0  I SCn ARGE  .  AND  TICAL 
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APPENDIX  D:  ERROR  MESSAGES  FOR  CE-QUAL-ELV2 


General 

CE-QUAL-ELV2  identifies  two  types  of  errors:  (1)  fatal  errors  which  stop 
the  computation,  (2)  non-fatal  errors  which  permit  the  computation  to 
continue.  The  messages  written  on  the  detection  of  these  errors  are 
discussed  below. 

1 .  "FATAL  INPUT  ERROR" 

This  error  occurs  when  the  input  stream  is  missing  a  required  card  or 
the  deck  is  out  of  order  (see  Appendix  A) .  The  code  checks  the  card 
identifiers  (field  0)  and  stops  the  computation  when  an  error  is 
detected.  An  input  error  also  occurs  when  too  few  or  too  many  tribu¬ 
taries  or  withdrawals  are  specified.  In  both  cases,  CE-QUAL-ELV2  writes  an 
error  message,  repeats  the  offending  card,  and  stops  the  computation. 

2.  "NON-FATAL  ERROR:  AZO  EXCEEDS  AZMAX" 
and 

"NON-FATAL  ERROR:  DZO  EXCEEDS  DZMAX" 

The  variables  AZO  and  DZO  are  set  by  the  user  and  represent  base  values 
of  vertical  dispersion  coefficients  for  the  x-momentum  and  for  tempera¬ 
ture,  salinity,  and  constituent  equations.  These  variables  are  used  with 
the  Richardson  number  to  compute  spatially  and  temporally  varying  disper¬ 
sion  coefficients.  The  upper  limits  on  the  latter,  however,  are  values 
based  on  the  horizontal  layer  thickness  h  and  the  time  step  At  (see 
Appendix  E) .  The  values  of  AZO  and  DZO  should  always  be  less  than 
these  maximum  values.  This  error  is  not  s  irious  enough  to  stop  the  com¬ 
putation,  however,  since  the  code  always  insures  that  the  computed  values 
of  AZ  and  DZ  do  not  exceed  their  computational  limits. 

3.  "FATAL  TVDS  ERROR  AT  _  DAYS" 

This  error  message  should  be  coded  in  the  TVDS  subroutine  to  detect  the 
case  where  CE-QUAL-ELV2  asks  for  time-varying  data  that  is  not  available-- 
for  example,  if  a  simulation  were  set  up  for  a  370-day  run,  with  only  one 
year  of  data  available.  The  illegal  time  request  is  printed  with  the 
error  message  and  the  computation  is  stopped. 


APPENDIX  E:  GLOSSARY  OF  IMPORTANT 
FORTRAN  VARIABLES 


The  glossary  lists  only  those  FORTRAN  variables  found  on  the  CE-QUAL-ELV2 
printed  output,  in  the  FORTRAN  statements  labelled  CHANGE,  in  the  TVDS  arg¬ 
uments,  and  in  this  report.  Those  variables  discussed  in  the  card  descrip¬ 
tions  (Appendix  A)  are  not  repeated  here. 

A:  vector  containing  subdiagonal  (backward)  coefficients  for  the  tri¬ 

diagonal  algorithm  solutions  of  Z,  T,  and  C 

n 

c  2  3  2 

ADMX:  ADvection  of  Momentum  in  the  X-direction:  y-  (U  Bh) ,  m  /s 

3  2 

ADMZ:  ADvection  of  Momentum  in  the  Z-direction:  (u,v7,  b),  ,  ,  m  /s 

b  b  k-% 

AR:  waterbody  area  by  layer,  m  (vector) 

AS:  A  Storage:  array  for  storing  A  vectors  for  use  in  constituent  con¬ 

centration  computations 

AX:  A^ ,  m  / s :  x-direction  momentum  dispersion  coefficient 

AZ :  Az,  m  /s:  dispersion  coefficient,  z-direction,  x-momentum  equation 

(array) 

AZMAX:  maximum  permitted  value  of  A z  (0.85*%-h‘ /At) ,  m2/s 

AZMIN:  minimum  value  of  A  (molecular  value),  m2/s  (scalar ) 

z 

AZO:  base  value  for  A  computation,  m2/s  (scalar' 

z 

B:  B,  m:  cell  width  (array) 

BETA:  (3,  fraction  of  incident  solar  radiation  absorbed  at  the  water 

surface 

BH1:  B'h,  m2 :  cross-sectional  area  of  a  cell  at  the  previous  time  step 

BH2:  B*h,  m2 :  cross-sectional  area  of  a  cell  at  the  current  time  step 

C:  vector  containing  superdiagonal  (forward)  coefficients  for  the  tri- 

diagonal  algorithm  solutions  of  Z  and  T 

CtiZY:  CHeZY  resistance  coefficient,  m  Ys 

CIN:  upstream  inflow  constituent  concentration,  mg  1  (vector) 


El 


CR:  Constituent  Right,  mg  1  1 :  right-hand  (open-water)  boundary  con¬ 
stituent  profile  (array) 

CS:  C  Storage:  array  for  storing  C  vectors  for  use  in  constituent  con¬ 
centration  computations 

CSHE:  k,  kinematic  coefficient  of  surface  heat  exchange,  mg  s-1 

CTRIB:  array  containing  tributary  constituent  concentrations,  mg  1 

CZ:  C* ,  wind  resistance  coefficient  (varies  with  wind  speed) 

Cl :  array  containing  constituent  concentrations  at  the  previous  time 
step,  mg  l-1 

C2:  array  containing  constituent  concentrations  at  the  current  time 
step,  mg  1_1 

D:  vector  containing  right-hand  side  for  the  tridiagonal  algorithm 
solution  of  Z  and  T 

1  3 

DG:  Density  Gradient:  p- -g^r  (PBh) ,  m3/s2  (This  gradient  term  in  the 
x-momentum  equation  includes  only  that  portion  of  hydrostatic 
pressure  due  to  density  differences  and  does  not  include  that  por¬ 
tion  due  to  changes  in  elevation.) 

£ 

DM:  Diffusion  of  Momentum:  (UBh) ,  m3/s2 

DX:  D  ,  m2/s:  dispersion  coefficient,  x-direction,  heat  and  constituent 

balance  equation  (array) 

DXI:  base  value  of  DX,  m2/s  (scalar) 

DZ:  D  ,  m2/s:  dispersion  coefficient,  z-direction,  heat  and  constituent 

balance  equation  (array) 

DZMAX:  maximum  permitted  value  of  D  (0.85*4'h2 /At) ,  m2/.;  (scalar) 

z 

DZMIN:  minimum  value  of  D  (molecular  value),  m2/s  (scalar) 

z 

DZO:  base  value  of  DZ  computation,  m2/s  (scalar) 

EL:  water  surface  elevation,  m 

ELTM:  ELapsed  TiMe,  days:  sum  of  NAt  and  TMSTRT,  time  since  beginning 

of  calendar  year 

ET:  E,  equilibrium  temperature  of  surface  heat  exchange,  C 

G:  g,  m/s2:  gravitational  constant 
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GAMMA:  y ,  m_  1  :  exponential  decay  constant  for  absorption  of  solar  radi¬ 
ation  with  depth 

HIN:  H  INitial,  m;  retains  original  value  of  h  (layer  thickness)  for 
computation  of  varying  h  in  top  layer  of  cells  from  Z  and  also 
used  in  reinitializing  a  layer  of  cells  immediately  after  a  layer 
has  been  added  or  subtracted  due  to  estuarine  volume  changes 

HN:  Hn,  °C  m3s  1  or  mg  1  1m3s  1:  net  rate  of  heat  or  constituent  addi¬ 
tion,  source  term  for  heat  balance  or  constituent  equation 

1  3 

HPG:  Horizontal  Pressure  Gradient:  p"  ^  (PBh) ,  m3/s2 

HI:  h,  m:  vertical  slice  thickness  at  the  current  or  previous  time 
step  (Only  the  top  layer  of  cells  has  varying  thickness  due  to 
elevation  changes;  however,  the  location  of  the  top  layer  of 
cells  may  vary  due  to  extreme  elevation  changes.) 

H2:  h,  m:  vertical  slice  thickness  at  the  current  or  previous  time 

step 

I:  counter  for  sweeping  across  (longitudinally)  the  estuary  grid  in 

x-direction  DO-loops 

IB:  I  Begin:  variable  locating  the  first  active  cell  in  each  line  for 

use  in  the  tridiagonal  solution  of  the  heat  balance  equation;  IB 
is  found  in  an  initial  program  step  and  stored  in  the  array  LC 

IE:  I  End:  variable  locating  the  last  active  cell  in  each  line  for 

use  in  the  tridiagonal  solution  of  the  heat  balance  equation;  IE 
is  found  in  an  initial  program  step  and  stored  in  the  array  LC 

IFLAG:  logic  flag  for  subroutine  TVDS;  if  IFLAG  =  0  ,  instruct  TVDS 

to  read  time-varying  boundary  condition  data;  if  IFLAG  =  1  , 
instruct  TVDS  to  provide  boundary  condition  data  for  a  partic¬ 
ular  time;  if  IFLAG  =  -1  ,  TVDS  error  and  computation  stopped 

IFORM:  control  variable;  if  IFORM  =  0  ,  terminal-oriented  display  is 
produced  (80-character  width);  if  IFORM  =  1  ,  full  width 
(132-character)  display  is  produced 

IL:  I  Left:  I-index  of  current  left-hand  segment 

IMAX:  I  MAXimum:  I-index  of  the  rightmost  segment 

IMAXM1:  I  MAXimum  Minus  1:  I-index  of  the  next- to- rightmost  segment 

IMAXM2:  I  MAXimum  Minus  2:  I-index  of  the  second -to- rightmost  segment 

ISC:  Integer  Segment  Coordinates:  vector  containing  K-index  of  bottom 

active  cell  of  each  segment 
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IY:  Integer  Year:  year  number 

JC:  constituent  counter 

K:  counter  for  sweeping  down  (vertically)  the  estuary  grid  in  z-direction 
DO-loops 

KB:  K  Bottom:  K-index  of  bottom  active  cell  in  a  particular  segment; 

KB  is  found  in  array  ISC 

KMAX:  K-MAXimum:  K-index  of  the  bottom  layer  of  cells 

KMAXMl:  K-MAXimum  Minus  1:  K-index  of  the  next  to  bottom  layer  of  cells 

KOUT:  computed  layer  number  at  which  each  tributary  enters  the  LAEM2 

grid  (vector) 

KT:  K-Top:  K-index  of  the  currently  active  top  layer  of  cells 

L:  integer  variable  denoting  active  cells  (L-l)  or  inactive  cells 

(L  =  0) 

LA:  X 

LC:  Layer  Coordinates:  array  containing  (1)  layer  number  (K) ,  (2)  IB, 
and  (3)  IE  for  each  active  layer  of  cells 

MNAC:  Maximum  Number  of  Active  Cells  (This  variable  may  be  used  to  more 

accurately  estimate  the  amount  of  computer  time  required  by  a  simu¬ 
lation,  see  paragraphs  109-111.) 

N:  counter  for  time  steps 

NLC:  Number  of  Layer  Coordinates  (This  variable  is  equal  to  the  number 
of  separate  layers  in  a  grid  and  may  be  used  to  set  the  first  di¬ 
mension  in  the  array  LC  .) 

P:  pressure,  N/m2  =  kg*m/s"2/m2 

P0:  P0 

PHI:  <t>,  degrees:  wind  direction 

QET :  Q  Evaporation  Total,  m3/s:  evaporation  rate 

3  -1 

QIN:  upstream  inflow  rate,  m  s 

QOUT:  vector  containing  downstream  outflow  rates 


QTRIB:  Q  TRIButary,  m3/s:  vector  containing  tributary  inflow  rates 
QWD:  Q  WithDrawal,  m3/s:  vector  containing  withdrawal  rates 
RHO:  p,  kg/m3:  density  of  water 

RHOA:  pa ,  kg/m3:  density  of  air  (used  in  wind  stress  computation  and 

assumed  constant  at  1.25  kg/m3) 

RR:  constituent  reaction  rate,  s  1 


S:  Stress:  (t  b).  ,  ra3/s 

z  k-Hj 

SAIN:  upstream  inflow  salinity,  ppt 


SR:  Salinity  Right,  ppt:  right-hand  (open-water)  boundary  temperature 
salinity  profile  (vector) 

SRO:  Solar  Radiation,  (°C’m3/s)/m2 


STRIB:  Salinity  of  TRIButary,  ppt:  vector  containing  tributary 

salinities 


SI:  S,  ppt:  array  containing  salinities  at  the  previous  time  step 

S2:  S,  ppt:  array  containing  salinities  at  the  current  time  step 

T:  vector  containing  the  results  of  the  layer-by-layer  solution  of  the 
implicit  temperature  equation 

TAPE5:  read  device 


TAPE6:  write  device 


TAPE61:  binary  write  device 

TD:  Td,  D:  dew  point  temperature,  °C 

TIN:  upstream  inflow  temperature,  °C 

TMSTRT:  time  of  year  in  Julian  days  when  simulation  starts,  read  in  as 

days,  converted  internally  to  seconds 

TR:  Temperature  Right,  °C:  right-hand  (open-water)  boundary  temperature 
profile  (vector) 

TTRIB:  Temperature  of  TRIButary,  °C:  vector  containing  tributary  temper¬ 
ature 
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Tl:  T,  °C:  array  containing  temperatures  at  the  previous  time  step 

T2;  T,  °C:  array  containing  temperatures  at  the  current  time  step 

U:  U,  m/s:  array  containing  x-direction  velocity  components 

-l 

UM:  cell  average  x-direction  velocity,  m/s 

V:  vector  containing  main  diagonal  (centered)  coefficients  for  the 
tridiagonal  algorithm  solution  of  Z  and  T 

VOL:  estuary  volume  based  on  initial  volume  and  sum  of  inflows  and 
outflows,  m3 

VOLE:  cumulative  waterbody  volume  by  layer,  m3  (vector) 

VOLP:  estuary  volume,  m3  (computed  from  geometry  and  water  surface 
elevations) 

VRATIO:  ratio  of  VOLP  to  VOL,  % 

VS:  V  Storage:  array  for  storing  V  vectors  for  use  in  constituents 
concentration  computations 

W:  ,  m/s:  array  containing  vertical  velocities 

WA:  W^,  m/s:  wind  speed 

WKT:  vector  containing  a  corrected  vertical  velocity  for  the  top  layer 
of  cells  that  permits  a  perfect  water  balance  in  that  layer 

WM:  cell  average  z-direction  velocity,  m/s 

ZR:  Z  Right,  m:  right-hand  (open-water)  boundary  water  surface  devia¬ 
tion 

Zl:  Z,  m:  vector  containing  deviation  of  water  surface  from  datum 

(top  of  active  layer  of  cells)  for  each  column,  positive  downwards 
(The  "1"  denotes  a  particular  time  level.) 

Z2:  Z,  m:  vector  containing  deviation  of  water  surface  from  datum  (top 

of  active  layer  of  cells)  for  each  column,  positive  downwards  (The 
"2"  denotes  a  particular  time  level.) 


APPENDIX  F:  USER  NOTES  FOR  PREPROCESSOR  GIN 


GIN  is  a  preprocessor  code  that  prepares  bathymetry  (BA)  cards  for  use 

by  CE-QUAL-ELV2 .  GIN  uses  GEDA  output  as  input  and  may  be  used  to  preview 

the  grid  geometry  prior  to  creating  the  BA  cards.  The  code  is  short  and  is 
documented  internally  with  comment  cards.  Following  are  the  most  impor¬ 
tant  features  of  GIN: 

a.  The  basic  grid  parameters  Ax  and  h  are  selected  during 
the  GEDA  run.  Output  from  the  final  GEDA  run  is  modified, 
eliminating  all  output  lines  preceding  "GEOMETRIC  MODEL 
FOR  UNSTEADY  FLOW  PROGRAMS . "  This  modified  GEDA  output  is 
the  only  input  to  GIN. 

b.  The  variable  IFORM  can  be  set  to  0  or  1  in  line  27  of  GIN. 

If  IFORM  =  0  ,  the  grid  geometry  is  previewed  in  a  form 

similar  to  that  presented  by  CE-QUAL-ELV2 .  If  IFORM  =  1  ,  the 
BA  cards  are  written.  The  user  will  generally  run  GIN  with 
IFORM  =  0  until  a  satisfactory  grid  is  established,  then 

with  IFORM  =1  as  a  final  step  before  working  with  CE-QUAL-ELV2 . 

c.  The  segment  order  may  be  reversed  by  changing  the  variable 
IMAGE  . 

d .  GIN  may  be  used  to  modify  a  grid  geometry  by  changing  some 
lines  of  the  code.  The  example  shows  some  lines  marked 
"CHANGE"  in  columns  73-80  in  DO-loop  170  that  modify  widths 
to  fit  better  the  given  elevation-area-volume  table.  The 
best  compromise  in  fitting  such  tables  may  be  to  match 
areas  in  the  surface  layers  and  volumes  elsewhere. 

e.  The  user  is  responsible  for  obtaining  a  legal  geometry  as 
outlined  in  paragraphs  27-29  of  this  report. 

f.  GIN  used  GEDA  output  in  feet  and  converts  all  dimensions  to 


APPENDIX  0:  POSTPROCESSOR  WE 


WE  is  a  postprocessor  code  that  uses  the  velocity  fields  generated  by 
CE-QUAL-KLV2  to  produce  vector  plots  like  those  shown  in  Figures  G1  and  G2. 
WE  uses  DISSPLA,  a  proprietary  software  package  created  by  Integrated 
Software  Systems  Corporation  of  San  Diego  and  available  on  many  time¬ 
sharing  systems.  WE  uses  a  file  created  by  CE-QUAL-ELV2  when  the  IPLOT 
parameter  on  the  PL  card  is  set  to  one.  The  code  is  short  and  docu¬ 
mented  internally  with  comment  cards.  However,  familiarity  with  DISSPLA 
subroutines  is  required  to  use  WE.  Following  are  some  of  the  most  im¬ 
portant  features  of  WE. 

a.  CE-QUAL-ELV2  writes  the  contents  of  the  TITLE  and  GE  cards  to 
TAPE61  in  binary  prior  to  writing  the  coordinates  of  all 

the  CE-QUAL-ELV2  cells  and  the  velocity  fields.  WE  reads 
this  information  on  TAPE51. 

b.  WE  computes  all  the  scaling  parameters  and  axis  labels 
from  information  written  by  CE-QUAL-ELV2 .  The  user  can 
change  the  axis  sizes  and  the  location  of  the  plot  on  the 
paper  by  changing  the  value  of  XAXIS  ,  ZAXIS  ,  XSTART  , 
and  ZSTART  (all  in  inches) . 

c.  The  velocity  components  are  translated  to  a  displacement 
from  the  cell  center  by  multiplying  each  component  by 
DLT  ,  a  time  in  seconds,  ^afferent  applications  yield 
different  velocity  magnitudes,  so  the  user  may  need  to 
change  the  value  of  DLT  .  It  is  converted  to  days  and 
identified  as  the  scaling  parameter  on  each  plot. 

d.  The  velocity  components  used  by  W  are  averaged  to  the 
cell  center  so  that  U(I,K)  =  (U(I  -  1,K)  +  U(I,K))/2 
and  W  =  (W(I,K  -  1)  +  W(I,K))/2  .  This  averaging  pro¬ 
cess  may  result  in  stepped  velocity  profiles  adjacent 
to  vertical  boundaries. 

e.  The  user  may  select  the  number  of  velocity  plots  by 
manipulating  the  PSTART  and  PEND  parameters. 


f.  WE  also  writes  all  the  plot  definition  variables  and 
title  information  to  TAPE6  as  a  debugging  aid. 

g.  The  previous  cycle  shown  in  Figure  62  is  created  by 
inserting  the  equation  for  ZR  in  TVDS  into  WE. 
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Figure  G2.  Water  surface  profile,  second  of  two  plots 
produced  by  WE  for  each  time  step 


APPENDIX  H:  JULIAN  DATE  CALENDAR 


in  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 
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